{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.2 Incompressible Navier-Stokes equations\n",
    "We are solving the unsteady Navier Stokes equations\n",
    "\n",
    "Find $(u,p):[0,T] \\to (H_{0,D}^1)^d \\times L^2$, s.t.\n",
    "\\begin{align}\n",
    "\\int_{\\Omega} \\partial_t u \\cdot v + \\int_{\\Omega} \\nu \\nabla u \\nabla v + u \\cdot \\nabla u v - \\int_{\\Omega} \\operatorname{div}(v) p &= \\int f v  && \\forall v \\in (H_{0,D}^1)^d, \\\\ \n",
    "- \\int_{\\Omega} \\operatorname{div}(u) q &= 0 && \\forall q \\in L^2, \\\\\n",
    "\\quad u(t=0) & = u_0\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import netgen.gui\n",
    "%gui tk\n",
    "import tkinter\n",
    "from math import pi\n",
    "from ngsolve import *\n",
    "from netgen.geom2d import SplineGeometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Schäfer-Turek benchmark\n",
    "We consider the benchmark setup of from http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html . The geometry is a 2D channel with a circular obstacle which is positioned (only slightly) off the center of the channel. The geometry:\n",
    "<center> ![title](http://www.featflow.de/media/dfg_bench2_2d/geometry.png) </center>\n",
    "\n",
    "The viscosity is set to $\\nu = 10^{-3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from netgen.geom2d import SplineGeometry\n",
    "geo = SplineGeometry()\n",
    "geo.AddRectangle( (0, 0), (2, 0.41), bcs = (\"wall\", \"outlet\", \"wall\", \"inlet\"))\n",
    "geo.AddCircle ( (0.2, 0.2), r=0.05, leftdomain=0, rightdomain=1, bc=\"cyl\")\n",
    "mesh = Mesh( geo.GenerateMesh(maxh=0.08))\n",
    "mesh.Curve(3)\n",
    "# viscosity\n",
    "nu = 0.001"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Taylor-Hood velocity-pressure pair\n",
    "We use a Taylor-Hood discretization of degree $k$. The finite element space for the (vectorial) velocity $\\mathbf{u}$ and the pressure $p$ are:\n",
    "\\begin{align}\n",
    " \\mathbf{u} \\in \\mathbf{V} &= \\{ v \\in H^1(\\Omega) | v|_T \\in \\mathcal{P}^k(T) \\}^2\\\\\n",
    " p \\in Q &= \\{ q \\in H^1(\\Omega) | q|_T \\in \\mathcal{P}^{k-1}(T) \\}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 3\n",
    "V = VectorH1(mesh,order=k, dirichlet=\"wall|cyl|inlet\")\n",
    "Q = H1(mesh,order=k-1)\n",
    "X = FESpace([V,Q])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The `VectorH1` is the product space of two `H1` spaces with convenience operators for identity, `grad` and `div`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Stokes problem for initial values\n",
    "Find $\\mathbf{u} \\in \\mathbf{V}$, $p \\in Q$ so that\n",
    "\\begin{align}\n",
    "\\int_{\\Omega} \\nu \\nabla \\mathbf{u} : \\nabla \\mathbf{v} - \\int_{\\Omega} \\operatorname{div}(\\mathbf{v}) p &= \\int \\mathbf{f}  \\cdot \\mathbf{v}  && \\forall \\mathbf{v} \\in \\mathbf{V}, \\\\ \n",
    "- \\int_{\\Omega} \\operatorname{div}(\\mathbf{u}) q &= 0 && \\forall q \\in Q,\n",
    "\\end{align}\n",
    "with boundary conditions:\n",
    " * the inflow boundary data ('inlet') with mean value $1$\n",
    " \\begin{align}\n",
    " u_x(0,y) &= 6 y(0.41-y)/(0.41)^2, \\quad \\int_{0}^{0.41} u_x(0,y) dy = 1 \\\\\n",
    " u_y(0,y) &= 0\n",
    " \\end{align}\n",
    " \n",
    " \n",
    " * \"do-nothing\" boundary conditions on the ouflow ('outlet') and\n",
    " * homogeneous Dirichlet conditions on all other boundaries ('wall')."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(X)\n",
    "velocity = CoefficientFunction (gfu.components[0])\n",
    "Draw(velocity,mesh,\"u\")\n",
    "Draw(gfu.components[1],mesh,\"p\")\n",
    "\n",
    "# parabolic inflow at bc=1:\n",
    "uin = CoefficientFunction((1.5*4*y*(0.41-y)/(0.41*0.41),0))\n",
    "gfu.components[0].Set(uin, definedon=mesh.Boundaries(\"inlet\"))\n",
    "\n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### inflow profile plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEWCAYAAACT7WsrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd81eX5//HXlR2SkEAIATIIhC0yw3JAUKs4sU5QVFxoreP707baVm1rh7W21WpdWHEiiBsnzjCUvWRDCGSwSVhJyL5+f+TYRkzIScg5n5NzrufjkcfjjPvzOW9uk3N5f8Z9i6pijDHG1CfI6QDGGGN8lxUJY4wxDbIiYYwxpkFWJIwxxjTIioQxxpgGWZEwxhjTICsSxi+IyDoRyXSzbW8RWSUiR0TkThF5SUT+5OGILUpE/iQi+0Vkt4ikikixiAS73ssSkZuczmj8Q4jTAYxpCap6UhOa/wr4WlUHAYjISx4J5SEikgrcA3RV1b2ul6MdjGT8mI0kTCDqCqxzOkRDvh8RHEcqUFinQBjjMVYkjF8Qke0icpbr8e9FZJaIvOI6pLRORDJc730FjAX+7TpE06uefd0sItkiUiQis0Wki+v1P4jIk67HoSJSIiKPup5HikiZiLSvZ3+ZIlIgIr9xHSLaLiJX13n/JRF5RkQ+FpESYKyIxLry7xORXBG5X0SCXP/Gz4EurvwviUiaiKiI1HtkQERuEJENInJAROaISNcT7G4TQKxIGH91ETATiANmA/8GUNUzgPnA7aoaraqb624kImcADwNXAJ2BXNd+AOYCma7Hw4DdwGjX81HAJlUtaiBPJ6ADkARcB0wVkd513r8K+DMQAywAngRige7AGOBa4HpV/QI4F9jpyj/5eJ0gIuOB3wCXAAmuf/uM421jTF1WJIy/WqCqH6tqNfAqMNDN7a4GpqnqClUtB34NjBKRNGAh0FNE4qktDi8ASSISTe0X+dxG9v2Aqpar6lzgI2oL0ffeV9VvVLUGqAQmAL9W1SOquh34B3CNm/+Gum4FHlbVDapaBfwFGGSjCeMuKxLGX+2u87gUiGjocMwxulA7egBAVYuBQiBJVY8Cy6gtCKOpLQrfAqfSeJE4oKoldZ7nuj7re/l1HncAQuvmcD1OciP/sboC/xKRgyJyECgCpJn7MgHIioQxP7ST2i9WAEQkCogHdrhemgucAQwGlrqenwMMB+YdZ7/tXPv6Xqrrs75Xdzrm/dSOJroe034HTZcP3KKqcXV+IlX122bsywQgKxLG/NAM4HoRGSQi4dQenlnsOuQDtUXhWmC9qlYAWcBNwDZV3dfIvv8gImEicjpwAfBmfY1ch8hmAX8WkRjXoaG7gdea8e95Fvi1iJwE4Dohfnkz9mMClN0nYUwdqvqFiDwAvA20o/Zw0oQ6Tb4FIvnfqGE9UMbxRxFQe/jrALWjh1LgVlXdeJz2d1B78jrHtf/ngWlN+scAqvqu65zJTFexOUTt1VH1FihjjiW26JAxnuW6E/w1VU12OosxTWWHm4wxxjTIioQxxpgG2eEmY4wxDbKRhDHGmAa1+qubOnTooGlpac3atqSkhKioqMYbGusrN1k/ucf6yT2e7Kfly5fvV9WExtq1+iKRlpbGsmXLmrVtVlYWmZmZLRvIT1lfucf6yT3WT+7xZD+JSG7jrexwkzHGmOOwImGMMaZBViSMMcY0yIqEMcaYBlmRMMYY0yArEsYYYxpkRcIYY0yDWv19EsY4RVUpLKmg4MBR8otKKThwlA1bKlhRsane9qHBQSS1iyS5XRtS2kfSMSaC4CDxcmpjmsaKhDFu2Hu4jIU5hazMO0heUel/i8LRyuoftBOAnOx693HsNGmhwUJS3P+KxslJcYxKjyctvg0iVjyMb7AiYUw99heXsyinkIVbC1mYU0jOvtrlqaPCgukaH0X3hChG90og5b8jgzYktYtk2cIFDd4hW1ZZzc6DR2tHHgdKfzAC+XTtbmYsqV3mulPbCEalxzOqezyj0uNJad/GW/9sY37EioQxLlv3FfPmsgK+2riHzXuKgdqiMKxbe67MSGFUejwndYlt9iGiiNBguidE0z0h+kfvqSo5+0v+W5Tmb9nHuytrl7ROiotkdK8ELhuaxJDUdjbKMF5lRcIEtCNllXz03S7eXF7A8twDBAcJo7rHM35QEqPS4zk5KZbQYM9f3yEipCdEk54QzaSRXVFVtuwtri0aWwt5f9UOZizJo3tCFJcPTeGSIUkkto3weC5jrEiYgFNToyzeVsSby/P5ZM1ujlZW06NjNL8+tw8/HZJExxjnv3xFhF6JMfRKjOG6U9IoLq/i4+928ebyfB75dCOPztlIZu+OXD40mTP7JhIWYhcqGs+wImECRmV1De+sKODprK3kFpYSEx7CxYOTuCIjmUEpcT59GCc6PIQrhqVwxbAUcvYV89byAt5eUcBXG/fSPiqM609JY/KpacREhDod1fgZKxLG71VV1/Deqp088eUW8opKGZAcy+NXDuKckzoRGRbsdLwm654Qza/G9eGes3szf8s+Xl2Yyz8+38wL32xjyujuXDcqjahw+9M2LcN+k4zfqq5RPli9k399uYVt+0s4qUtbXrgugzP6dPTpUYO7goOEzN4dyezdke8KDvLY55v526eb+M/8bdwyujvXjOpKmzD7Ezcnxn6DjN+pqVE+WrOLx7/YzNZ9JfTpFMNz1wzl7H6JflEc6jMgOY4Xrx/OyrwDPPbFFh7+ZCPPz8/h1jHpTBrZlYjQ1jdiMr7BioTxKxt2Hea+d9awOv8gvRKjeebqIZxzUieCAuTO5sGp7XjlhuEs217EY19s5k8fbeDFb7bzp5/2Z2zvjk7HM62Q1y6JEJFpIrJXRNY20m6YiFSJyGXeymZav7LKah75dCMXPrmAgqJS/nnFQD69azTnntw5YApEXRlp7Zl+00hm3DySyLBgrn9xKXfMWMm+I+VORzOtjDevm3sJGHe8BiISDDwCfOaNQMY/fJO9n3Men8czWVv56eAkvrxnDJcMSQ7I4nCsUenxfHTnafy/s3oxZ+1uzvrnXGYtzUePnSPEmAZ4rUio6jygqJFmdwBvA3s9n8i0dgdKKrhn1mqu/s9iBHj9phE8evlA4tqEOR3Np4SHBHPXWT35+K7T6J0Yw6/e/o6Jzy8iZ1+x09FMKyDe/D8KEUkDPlTV/vW8lwS8DowFprnavdXAfqYAUwASExOHzpw5s1l5iouLiY7+8RQJ5sd8ra8W7axi+oZySqvg3G6hXJQeSliw8yMHX+unY9WoMq+gijc2VVBZAxelh3J+t1Cvz0br6/3kKzzZT2PHjl2uqhmNtfOlE9ePA/eqak1jV6Co6lRgKkBGRoY2NKFaY7KyshqcjM38kK/01dGKah58fy1vflfAoJQ4/nrpyfTp1NbpWP/lK/10PGcAtx0u4w8frOedNbsoqIrh3xMH09GL03y0hn7yBb7QT75UJDKAma4C0QE4T0SqVPU9Z2MZX5Gzr5jbpq9g054j3HlGD+46q5etx9BMHdtG8NTVQzhzRQG/fXct5z0xnycmDOaUHh2cjmZ8jM9M+KKq3VQ1TVXTgLeA26xAmO99+N1OLnxyAXsOl/Hi5GHcfXZvKxAt4JIhybx/+6nERoYy6YXFPPnlFmpq7KS2+R+vjSREZAaQCXQQkQLgd0AogKo+660cpnUpr6rmLx9t4OWFuQxJjePfVw2hS1yk07H8Sq/EGGbffhq/eXcN//h8M0tzD/D4lYNoH2UXABgvFglVndiEtpM9GMW0EvlFpdz++gpWFxzixtO6cd+5fbwybXcgigoP4fErBzEsrT0PfbCe85+Yz7+vGsLQru2cjmYcZn9xxict2LKfC55cQM6+Ep6dNJQHLuhnBcLDRIRJI7vyzm2nEBIsXPncQl5duN3pWMZh9ldnfM47KwqY/OISOrWN4MM7T2Nc/05ORwoo/ZNi+fCO0xnTK4EH3l/HXz/ZaOcpApgVCeMzVJWns7K5e9ZqhqW1582fjaJrfJTTsQJSbGQoz10zlKtHpPLs3K3cPWsVFVU1TscyDvClS2BNAKuuUX4/ex2vLsrlooFdePTyAYSH2MylTgoJDuJPF/enS1wkj87ZxL7icp6dNNQWNgowNpIwjiurrOZnry3n1UW53DK6O49fOcgKhI8QEX4+tgd/v3wgi3OKuPzZhew5XOZ0LONFViSMow6UVHDV84v4fMMefndhP359Xl+bmM8HXTY0mWmTh5FfVMolT3/Llj1HnI5kvMSKhHFMflEplz7zLWt3Hubpq4Zw/andnI5kjmN0rwTeuGUU5VU1XPrMtyzZ1th8ncYfWJEwjti2v4TLnv2WwpIKXrtxBOee3NnpSMYN/ZNiefe2U+gQE86kFxYzf8s+pyMZD7MiYbwut7CEiVMXUVmtvHHLSIZ3a+90JNMEKe3b8Natp9C9QxQ3vbyMb7P3Ox3JeJAVCeNV+UWlTJy6iLKqal67cYRPzeBq3Nc+KozpN42ga3wbbnh5KYtyCp2OZDzEioTxmoIDpUyYuoiSitoC0a+LFYjWLD46nOk3jSS5XRtueGmpnaPwU1YkjFfsPHiUic8v4khZJa/dOIL+SbFORzItICEmnNdvHkGn2Aiuf3EJy3OtUPgbKxLG43YfKmPi84s4WFLJqzeO4ORkKxD+pGNMBDNuHknHthFcN20pK/MOOB3JtCArEsaj9hyuLRCFxRW8fONwBqbEOR3JeEBi29pCER8dxrUvLGF1/kGnI5kWYkXCeMz+4nImPr+IvYfLePmGYQxJtWmn/Vmn2NpCERcVyjUvLGbtjkNORzItwIqE8YjSiipufGkpOw8e5cXrhzO0q13mGgi6xEUy4+aRRIeHcP1LSyk4UOp0JHOCrEiYFlddo9w5YxVrdhziiQmD7T6IAJPcrg0v3TCcsspqJr+4lEOllU5HMifAioRpUarKHz5Yxxcb9vD7i07i7JNsLYhA1CsxhqnXZJBbWMKUV5dRXlXtdCTTTF4rEiIyTUT2isjaBt6/WkS+E5E1IvKtiAz0VjbTcp6fn8MrC3OZMro7145KczqOcdCo9Pja2WO3FfHLN7+zhYtaKW+OJF4Cxh3n/W3AGFU9GfgjMNUboUzL+fC7nfzl442cP6Az943r43Qc4wPGD0ril+f0Zvbqnfz9s01OxzHN4LVFh1R1noikHef9b+s8XQQkezqTaTlLthVx9xurGZbWjn9cPtCm+zb/dVtmOjsOHuXprK0ktYvk6hFdnY5kmkBUvTcEdBWJD1W1fyPtfgH0UdWbGnh/CjAFIDExcejMmTOblae4uJjo6OhmbRtojtdXO4tr+PPio8SECfePiCQ6LHALhP1O1a+6RvnXynLW7KvmriHh9GhTZv3kBk/+Po0dO3a5qmY01s7nioSIjAWeBk5T1UZnDcvIyNBly5Y1K09WVhaZmZnN2jbQNNRX+46U89Onv6Gsspp3bzuVlPZtvB/Oh9jvVMNKyquYMHUR2XuL+VVGKNePP9PpSD7Pk79PIuJWkfCpq5tEZADwH2C8OwXCOKu8qpopry6jsLiCaZOHBXyBMMcXFR7CC5MzaB8VxuMrytlry6C2Cj5TJEQkFXgHuEZVNzudxzTuoQ/WszLvIP+8YiADkm26DdO4jjERTJs8jKNVym3TV1BRVeN0JNMIb14COwNYCPQWkQIRuVFEbhWRW11NHgTigadFZJWINO8YkvGKWUvzmb44j1vHpNuqcqZJeneK4cb+4SzLPcCfP1rvdBzTCG9e3TSxkfdvAuo9UW18y3cFB7n//bWc1qMDvzi7l9NxTCs0onMIVW278Pz8bQxIjuPSoXYxo6/ymcNNpnUoKqngZ6+tICE6nCcmDiYk2H6FTPPcO64PI7u35zfvrrHJAH2Y/YUbt1VV13DHjBXsKy7nmUlDaB8V5nQk04qFBAfx76uG0K5NGLe+tpyDpRVORzL1sCJh3Pb3zzbzTXYhf7q4v52oNi2iQ3Q4z0wawt7D5dw5cxXVNnWHz7EiYdyydHcVz87dytUjUrkiI8XpOMaPDE5txx/Gn8S8zft47HO7sNHXWJEwjdqy5wgvrClncGocD17Yz+k4xg9NHJ7KhGEp/PvrbOas2+10HFOHFQlzXMXlVdzy2nLCguGZq4cSHhLsdCTjp35/0UkMTI7lnlmrydlX7HQc42JFwhzX72evY/v+En42MIJOsRFOxzF+LCI0mGcmDSUkWLhr5iq70c5HWJEwDfpg9U7eWl7Az8f2oG+8jSCM53WJi+SvlwxgzY5D/NPOT/gEKxKmXgUHSvnNu2sYnBrHnWf2dDqOCSDj+ndi4vBUnpu3lW+y9zsdJ+BZkTA/UlVdw/97YxWq8K8rBxNqN8wZL3vggr506xDF3bNWcaDE7p9wkv31mx95OmsrS7cf4I8Xn0RqvM3saryvTVgIT0wYTFFJBfe+/R3eXNLA/JAVCfMDy3MP8K8vt3DxoC78dLDNp2Oc0z8plnvH9eGz9Xt4fUme03EClhUJ81+Hyyq5a+ZKOsdG8NDFx1080BivuOHUbpzeswN//HA92XuPOB0nIFmRMP/14Htr2XWojH9NGEzbiFCn4xhDUJDwj8sH0iYshDtmrKK8qtrpSAHHioQB4N2VBby3aid3ntGToV3bOR3HmP/q2DaCRy8bwIZdh3n0001Oxwk4ViQMeYWlPPDeOoaltePnY9OdjmPMj5zZN5FrR3XlPwu2MXfzPqfjBBQrEgGupkb5xVurEeCxKwfZ+hDGZ/3mvL70SozmV2+t5tDRSqfjBAxvLl86TUT2isjaBt4XEXlCRLJF5DsRGeKtbIFs+uJclmwr4oEL+pHczi53Nb4rIjSYv18+kP3FFfzlow1OxwkY3vzfxpeAccd5/1ygp+tnCvCMFzIFtPyiUh7+ZCOn9+zA5Rl2uavxfQOS45gyujtvLMtn/hY77OQNXisSqjoPKDpOk/HAK1prERAnIp29ky7wqCq/eXcNAjx8ycmIiNORjHHLXWf2pHtCFPe9vYbi8iqn4/g98eadjCKSBnyoqj+6CF9EPgT+qqoLXM+/BO5V1WX1tJ1C7WiDxMTEoTNnzmxWnuLiYqKjo5u1bWs3t6CSF9dWcG2/MM5Ibfxy10Duq6awfnLPifbTlgPV/GVxGWNTQ7i2X3gLJvMtnvx9Gjt27HJVzWisXYhHPt3DVHUqMBUgIyNDMzMzm7WfrKwsmrtta7b7UBl3fD2XEd3a8/tJIwkKanwUEah91VTWT+450X7KBHaFrGfaN9uYMi6Dkd3jWyqaT/GF3ydfupRlB1B3Xcxk12umBakqv313DZU1NTxy6QC3CoQxvugX5/QitX0b7n37O45W2E12nuJLRWI2cK3rKqeRwCFV3eV0KH/z/qqdfLlxL784uzdpHaKcjmNMs7UJC+GRSweQW1jKPz6zm+w8xWuHm0RkBrWjxA4iUgD8DggFUNVngY+B84BsoBS43lvZAsW+I+X8/oN1DEmN4/pTuzkdx5gTNio9nkkjU3nhm22cN6AzQ1JttoCW5rUioaoTG3lfgZ97KU5AevD9tZRWVPO3ywYSbIeZjJ+479y+fL1xH798czUf3Xk6EaG2imJL8qXDTcaDPl6zi0/W7ub/zupJj4529Y3xH9HhIfzlkpPZuq+EJ77c4nQcv2NFIgAcKq3kwffXcnJSLFNO7+50HGNa3JheCVw+NJnn5uWwbuchp+P4FSsSAeDvn22iqKSChy852eZmMn7r/vP7ERsZygPvraWmxlayayn2jeHn1hQc4rXFuVw7Ko3+SbFOxzHGY2LbhPLrc/uwIu8gby0vcDqO37Ai4ceqa5T731tDfFQ4d5/dy+k4xnjcpUOSyejajoc/2cCBkgqn4/gFKxJ+bObSPFYXHOK35/exleZMQAgKEv54cX8Ol1Xxtzl270RLsCLhpwqLy/nbp5sY0a09Fw9KcjqOMV7Tt3NbJp+SxsyleazKP+h0nFbPioSfeuTTjZSUV/Gni/vbDK8m4PzfWT1JiA7n/vfWUG0nsU+IFQk/tGx7EbOWFXDjad3omRjjdBxjvC4mIpT7L+jH2h2Hmb441+k4rZoVCT9TVV3D/e+tpXNsBHee2dPpOMY45sIBnTm1RzyPztnEviPlTsdptaxI+JlXFuaycfcRHrygH1HhrXImeGNahIjw0Pj+lFVW8/AnttxpczW5SIhIlIjY5Cg+aM/hMv75+WZG90pgXP9OTscxxnHpCdHcfHp33lmxg8U5hU7HaZUaLRIiEiQiV4nIRyKyF9gI7BKR9SLyqIj08HxM444/f7SBiuoaHrroJDtZbYzLHWf0JCkukgfeX0tldY3TcVodd0YSXwPpwK+BTqqaoqodgdOARcAjIjLJgxmNGxZuLWT26p3cOibd1okwpo7IsGB+d2E/Nu8p5uVvtzsdp9Vx56D1WapaeeyLqloEvA28LSJ2p5aDqmuUP364nqS4SG7LTHc6jjE+5yf9EhnTK4F/fbmFS4Yk0z4qzOlIrUajI4nvC4SIfC4iA4/Xxjjj7eUFrN91mHvP7WNz6RtTDxHh/vP7UlpRzeNfbHY6TqvSlBPX9wKPi8iLItLZU4FM05SUV/HoZ5sYnBrHhQPsP4sxDemZGMPE4SlMX5xH9t4jTsdpNdwuEqq6QlXHAh8Cn4rI70QksikfJiLjRGSTiGSLyH31vJ8qIl+LyEoR+U5EzmvK/gPRc3O3su9IOQ9c0M9OVhvTiP93Vi/ahAbzl483Oh2l1WjSJbBS+y20CXgGuAPYIiLXuLltMPAUcC7QD5goIv2OaXY/MEtVBwMTgKebki/Q7Dx4lKnzc7hwYBdb29cYN8RHh3P7GT34auNe5m/Z53ScVsHtIiEi3wA7gMeAJGAykAkMF5GpbuxiOJCtqjmqWgHMBMYf00aBtq7HscBOd/MFor/P2USNwr3jejsdxZhWY/KpaaS0j+TPH22weZ3cIKrudZKInASs13o2EJENqtq3ke0vA8ap6k2u59cAI1T19jptOgOfAe2AKGqvrFpez76mAFMAEhMTh86cOdOtf8OxiouLiY5unes9bztUzR8WlnFB91Au6+X5KzVac195k/WTe5zup6W7q3hqVTnXnxTGmBTfvTjTk/00duzY5aqa0Vg7t+dtUNV1x3n7fHf304iJwEuq+g8RGQW8KiL9VfUHd8Co6lRgKkBGRoZmZmY268OysrJo7rZOUlWeem4hHaJrePjaTGK8sFZEa+0rb7N+co/T/TRGlcUHFvJBbin3XHEa0T46hY3T/QQtNHeTqua40WwHkFLnebLrtbpuBGa59rkQiAA6tERGf/Lp2t0s3X6Au3/S2ysFwhh/IyLcf0E/9heX80xWttNxfFqzi4SIdBaR8CZsshToKSLdRCSM2hPTs49pkwec6dp/X2qLhJ1dqqO8qpqHP9lI78QYrshIdjqOMa3WoJQ4Lh7Uhefnb6PgQKnTcXzWiYwkXgU2isjf3WmsqlXA7cAcYAO1VzGtE5GHROQiV7N7gJtFZDUwA5hc3zmQQPbKt7nkFZVy/wV9CQm2SXyNORG/GtcHAR61pU4b1OwDcap6luuS2GMvYz3eNh8DHx/z2oN1Hq8HTm1uJn9XVFLBE19tYWzvBE7vmeB0HGNavS5xkUwZ3Z0nv8pm8ilpDLZLyX/khP5XVGsd74S2aUH/+mIzpRXV/Oa8415IZoxpglvHpJMQE86fPtqAHbj4MbeKhIj0EZF7ReQJ18+9rnMGxkvyCkuZvjiPCcNSbElSY1pQVHgI9/ykF8tzD/DFhr1Ox/E57qwncS+1N74JsMT1I8CM+qbWMJ7x2BebCQkW7rIlSY1pcZcNTaZ7hyj+PmeT3WB3DHdGEjcCw1T1r6r6muvnr9TeQX2jZ+MZgI27D/Peqh1MPqUbHdtGOB3HGL8TEhzE3Wf3YtOeI8xefeyV+YHNnSJRA3Sp5/XOrveMh/3js81Eh4dw65juTkcxxm+d178z/Tq35bHPt1BRZV9t33OnSPwf8KWIfCIiU10/nwJfAnd5Np5ZkXeAz9fv4ZbR3YlrYwulGOMpQUHCL8f1Jq+olDeW5Tsdx2c0egmsqn4qIr2oPbyU5Hp5B7BUVas9GS7QqSqPfrqJDtFhXH9qN6fjGOP3MnslMDytPU9+uYXLhiQTGWaLeLlz4lpUtUZVF6nq266fRXULhNhCBh7xTXYhC3MK+fnYHkT56NwyxvgTkdrRxN4j5by8cLvTcXyCO4ebvhaRO0Qkte6LIhImImeIyMvAdZ6JF7hUlUfnbCQpLpKrRqQ2voExpkUMS2vP2N4JPJO1lUNHbWVmd4rEOKCa2kted4rIehHZBmyhdtbWx1X1JQ9mDEhz1u1hdcEh7jqrJ+EhNuQ1xpvuObs3h45W8p/57sxd6t/cOSdRRu0KcU+LSCi1s7JWqup+T4cLVNU1yj8+20R6QhSXDE5qfANjTIvqnxTLBQM688KCbVw7Ko2EmKbMZepfmjQth6pWquou4K/fr28tIqM9kiyAvbtyB1v2FvOLs3vbJH7GOOTun/SivKqGp74O7KnEm/sN9DvgRRF5FRjWgnkCXnlVNY99vpmTk2IZ17+T03GMCVjdE6K5fGgyry/OC+ipxJtbJP4IbKR2TepZLRfHzFySz46DR/nlOb2xi8aMcdadZ/YEgX99scXpKI5pbpH4lar+HriN2lGFaQGlFVU8+VU2I7q15/SetiCfMU7rEhfJNSO78vaKArL3HnE6jiPcLhIiskxEXhCR/wMGikiCqhYDt3guXmCZviiP/cXl/GqcjSKM8RW3ZaYTERrMk18F5rmJpowkLgLeBMKoLQzbRSTX7rpuGUcrqnlu3lZO79mBoV3bOx3HGOMSHx3OtaPS+GD1TrbuK3Y6jte5XSRUdaeqfqqqf1PVK4AM4D+eixZYXl+Sx/7iitpjoMYYn3LT6d0IDwnmqQAcTTTlcFPXus9VdQPQqykfJiLjRGSTiGQ3tBaFiFzhumFvnYi83pT9t1ZlldU8O3cro7rHMyzNRhHG+JoO0eFMGpnKe6t2sH1/idNxvKoph5tmiEiBiMwXkadF5C9Af3c3FpFg4CngXGrXxZ4oIv0zVVpmAAATiklEQVSOadMT+DVwqqqeRO0MtH7vjaX57DtSbqMIY3zYzaO7ExocFHD3TTTlcNMpQApwPfA5sBW4sAmfNRzIVtUcVa2gdrW78ce0uRl4SlUPuD7T79cSLK+q5pmsrQxPa8/I7jaKMMZXdYyJ4KoRqbyzcgd5hYFz30STphbV2lXCs10/TZUE1J2kvQAYcUybXgAi8g0QDPxeVT89dkciMgWYApCYmEhWVlYz4kBxcXGzt20pX+VVsvtwBdf2hrlz5zqa5Xh8oa9aA+sn97TWfhoQWoOg3D9jPjf09/xUHb7QT742/3QI0BPIBJKBeSJysqoerNtIVacCUwEyMjI0MzOzWR+WlZVFc7dtCRVVNfz271kM7dqOn10yyqcve3W6r1oL6yf3tOZ+WlW+lumL8/jL1cNJbtfGo5/lC/3kzYmBdlB7uOp7ya7X6ioAZrvmiNoGbKa2aPilt1cUsOPgUe48s6dPFwhjzP/cmplOkAjPZG11OopXeLNILAV6ikg3EQkDJgCzj2nzHrWjCESkA7WHn/xyrt7K6tqJwwamxDHa7q42ptXoHBvJ5RnJzFqWz86DR52O43FeKxKqWgXcDswBNgCzVHWdiDwkIhe5ms0BCkVkPfA18EtVLfRWRm96d+UOCg4c5a4ze9gowphW5meZ6ajCs3P9fzTh1XMSqvox8PExrz1Y57ECd7t+/FaVaxRxclIsY3t3dDqOMaaJktu14bKhycxcks9tmT3oFBvhdCSPscUKHDB79U5yC0vtXIQxrdjPx/agWpXn5vn3aMKKhJdV1yj//iqbvp3bclZfG0UY01qltG/DJYOTeH1xHnuPlDkdx2OsSHjZh9/tJGd/iZ2LMMYP/HxsDyqra5g61y+vrwGsSHiVqvL011vplRjN2f1s1TljWru0DlGMH5TE9MV5HCipcDqOR1iR8KKvN+1l054j3DomnaAgG0UY4w9uHZPO0cpqXl2U63QUj7Ai4UXPzs0hKS6SCwd2cTqKMaaF9O4Uwxl9OvLSt9s5WuF/y+tYkfCSFXkHWLKtiBtP60ZosHW7Mf7k1jHpFJVU8Oby/MYbtzL2beUlz2ZtJTYylCuHpTTe2BjTqgxLa8fg1DimzsuhqrrG6TgtyoqEF2TvLebzDXu4blRXosJ9bU5FY8yJEhFuHZNOwYGjfLx2t9NxWpQVCS+YOm8r4SFBXHdKmtNRjDEe8pO+iaQnRPFs1lZqJ4/wD1YkPGzP4TLeXbmDKzJSiI/2/PzzxhhnBAUJt4xOZ/2uw8zfst/pOC3GioSHTVuwjeoa5ebTuzsdxRjjYeMHdyGxbbhfTfxnRcKDDh2tZPriPM4f0IWU9p5dnMQY47zwkGBuOLUb324t5LuCg41v0ApYkfCg6YtzKS6v4pbRNoowJlBcNSKVmIgQnvOTqTqsSHhIWWU10xZs5/SeHeifFOt0HGOMl8REhDJpZFc+WbuL7ftLnI5zwqxIeMg7K3awv7icn41JdzqKMcbLrj81jZCgIKbOb/2jCSsSHlBdozw/P4eTk2IZlR7vdBxjjJd1jIng0qFJvLW8oNVPI+7VIiEi40Rkk4hki8h9x2l3qYioiGR4M19L+WzdbrbtL+HWMek2HbgxAerm07tTWV3Dy99udzrKCfFakRCRYOAp4FygHzBRRPrV0y4GuAtY7K1sLUlVeXbuVtLi2zCuv00Hbkyg6p4QzbiTOvHqwlyOlFU6HafZvDmSGA5kq2qOqlYAM4Hx9bT7I/AI0CrHaEu3H2B1wSFuOr07wTYduDEB7ZYx6Rwuq+LNZQVOR2k2b04klATUnSKxABhRt4GIDAFSVPUjEfllQzsSkSnAFIDExESysrKaFai4uLjZ2zbkyZVlRIVCQkkOWVnbWnTfTvJEX/kj6yf3BFI/9YwL4pkvN5BWuZ2gJh5+9oV+8pnZ5kQkCPgnMLmxtqo6FZgKkJGRoZmZmc36zKysLJq7bX3yi0pZOedrbh2Tzjln9mmx/fqClu4rf2X95J5A6qfS+F3cNn0FVR37cvZJTTsE7Qv95M3DTTuAuvNkJ7te+14M0B/IEpHtwEhgdms6ef3yt7X/p3DNqK5ORzHG+Iiz+yWSFBfJtG9a55EFbxaJpUBPEekmImHABGD292+q6iFV7aCqaaqaBiwCLlLVZV7M2GzF5VW8sTSf807uTOfYSKfjGGN8REhwENed0pVFOUWs23nI6ThN5rUioapVwO3AHGADMEtV14nIQyJykbdyeMqby/I5Ul7FDad1czqKMcbHXDkslTZhwUxbsN3pKE3m1XMSqvox8PExrz3YQNtMb2RqCdU1ykvfbmdIahyDUuKcjmOM8TGxkaFcPjSZGUvyuffc3nSMiXA6ktvsjusW8NXGveQWltoowhjToMmndqOiuobpi/KcjtIkViRawLQF2+gSG8G4Jl65YIwJHN06RHFmn45MX5xLWWW103HcZkXiBK3feZiFOYVcd0oaIcHWncaYht1wWjf2F1fwweqdTkdxm32rnaBp32wjMjSYCcNSnY5ijPFxp6TH0zsxhhcWbGs162BbkTgB+46UM3vVTi4bmkxsm1Cn4xhjfJyIcMNpaWzcfYSFOYVOx3GLFYkTMH1xLhXVNUw+Nc3pKMaYVmL8oCTaR4W1msthrUg0U3lVNa8tymVs7wTSE6KdjmOMaSUiQoOZNCKVLzfuaRUr11mRaKYPVu9if3GFXfZqjGmySSO7EhIkvNQK1pqwItEMqsq0Bdvo2TGa03p0cDqOMaaV6dg2ggsHdOHNZfkc9vG1JqxINMOinCLW7zrMDad1s5XnjDHNcv2p3SipqGbW0vzGGzvIikQzvLJwO3FtQvnp4CSnoxhjWqmTk2MZntaeVxbmUlPju5fDWpFool2HjvLZ+j1cmZFCRGiw03GMMa3YNaO6kldUytwt+5yO0iArEk00Y0k+NapcNcJunjPGnJhzTupEh+hwXluY63SUBlmRaILK6hpmLMljTK8EusZHOR3HGNPKhYUEMXF4Cl9t2kt+UanTceplRaIJPlu3h31HyrnWVp4zxrSQicNTEeD1Jb45O6wViSZ4ZeF2kttFMqZXR6ejGGP8RJe4SH7SL5E3lub75OywViTctHnPERZvK+LqEV0JDrLLXo0xLeeakWkUlVTwydpdTkf5Ea8WCREZJyKbRCRbRO6r5/27RWS9iHwnIl+KiM8c13ltUS5hIUFcOSzF6SjGGD9zSno83TtE8aoPnsD2WpEQkWDgKeBcoB8wUUT6HdNsJZChqgOAt4C/eSvf8RSXV/HOih1ccHJn2keFOR3HGONngoKEq0d2ZUXeQdbuOOR0nB/w5khiOJCtqjmqWgHMBMbXbaCqX6vq96f4FwHJXszXoHdX7qC4vIpJdsLaGOMhlw1NJiI0iNcW+dZowptFIgmoe/95geu1htwIfOLRRG5QVV5bmEv/pLYMTolzOo4xxk/FRoZy8aAk3lu1g0NHfWc+pxCnA9RHRCYBGcCYBt6fAkwBSExMJCsrq1mfU1xc3Oi2m4qq2bSnjOv7hzF37txmfY4/cKevjPWTu6yf6tc3tJqZlTU88kYWZ6eF+kQ/ebNI7ADqnvVNdr32AyJyFvBbYIyqlte3I1WdCkwFyMjI0MzMzGYFysrKorFt33p9BTER+/jVFWcQGRa403C401fG+sld1k8Ne7/gGxbtr+RP145h3ry5jveTNw83LQV6ikg3EQkDJgCz6zYQkcHAc8BFqrrXi9nqtfdIGZ+u3c3lQ1MCukAYY7znmlFdydlfwrdbfWN5U68VCVWtAm4H5gAbgFmquk5EHhKRi1zNHgWigTdFZJWIzG5gd17xxpJ8qmqUSSNtniZjjHec27/2KspXF213Ogrg5XMSqvox8PExrz1Y5/FZ3sxzPFXVNby+JI/Te3aguy1PaozxkojQYK7ISGHqvK2ckxDpdBy747ohX2zYy65DZUwaaZe9GmO86+oRqSjwdX6V01GsSDTktUW5dImN4Mw+Nk+TMca7Utq34YzeHZmbX0VFVY2jWaxI1CO3sIQF2fuZMDyVkGDrImOM900a2ZXDFcoXG/Y4msO+Aesxc2k+QQJXZNg8TcYYZ4zulUD7CGGGw1OIW5E4RmV1DW8uK+CMPh3pFBvhdBxjTIAKDhJOTwphQfZ+RxcksiJxjK827mV/cTkThtllr8YYZ41Orr0A9c1l+Y209BwrEseYuSSPxLbhZPZOcDqKMSbAxUcGMaZXArOWFVBV7cwJbCsSdew8eJS5m/dxRUaKnbA2xviECcNS2X24jLmb9zny+fZNWMesZfkodsLaGOM7zuzbkQ7R4cxY4swhJysSLtU1yqyl+ZzWowMp7ds4HccYYwAIDQ7i8oxkvt60lz2Hy7z++VYkXOZt2cfOQ2V2wtoY43OuzEihukYdOYFtRcJl5pI84qPC+Em/RKejGGPMD6R1iGJU93jeWJZPTY169bOtSFA7JfiXG/Zy6dBkwkKsS4wxvmfC8BTyi47yzdb9Xv1c+0YE3lpeQFWNcuUwO2FtjPFN55zUibg2ocz08gnsgC8SNTXKG0vzGd6tPek2JbgxxkdFhAZzyeBkPlu/m8Liehft9IiALxKLcgrJLSxl4nAbRRhjfNvE4SlUVitvryjw2mcGfJGYsTSfthEhnNu/s9NRjDHmuHomxjC0aztmLs1H1TsnsAO6SBRXKHPW7uaSIclEhNoa1sYY3zdhWAo5+0pYuv2AVz7Pq0VCRMaJyCYRyRaR++p5P1xE3nC9v1hE0jyZ55udVVRU1zDBDjUZY1qJ8wd0JiY8hJlemkLca0VCRIKBp4BzgX7ARBHpd0yzG4EDqtoDeAx4xFN5VJW5BZUMSomjT6e2nvoYY4xpUW3CQhg/uAsfrdnFodJKj3+eN0cSw4FsVc1R1QpgJjD+mDbjgZddj98CzhQR8USYFXkH2FmsdsLaGNPqTBiWSnlVDe+t2uHxzwrx+Cf8TxJQ9wLfAmBEQ21UtUpEDgHxwA/uHhGRKcAUgMTERLKyspocZsuBavrGKbGHtpKVldPk7QNNcXFxs/o50Fg/ucf6yT3H66eRnYPZtT2brIrtHs3gzSLRYlR1KjAVICMjQzMzM5u8j0ygZ1YWzdk2EGVZX7nF+sk91k/uOV4/eav7vHm4aQdQ99hOsuu1etuISAgQCxR6JZ0xxpgf8WaRWAr0FJFuIhIGTABmH9NmNnCd6/FlwFfqrYuBjTHG/IjXDje5zjHcDswBgoFpqrpORB4ClqnqbOAF4FURyQaKqC0kxhhjHOLVcxKq+jHw8TGvPVjncRlwuTczGWOMaVhA33FtjDHm+KxIGGOMaZAVCWOMMQ2yImGMMaZB0tqvMBWRfUBuMzfvwDF3c5sGWV+5x/rJPdZP7vFkP3VV1YTGGrX6InEiRGSZqmY4naM1sL5yj/WTe6yf3OML/WSHm4wxxjTIioQxxpgGBXqRmOp0gFbE+so91k/usX5yj+P9FNDnJIwxxhxfoI8kjDHGHIcVCWOMMQ0KiCIhIuNEZJOIZIvIffW8Hy4ib7jeXywiad5P6Tw3+mm0iKwQkSoRucyJjL7Cjb66W0TWi8h3IvKliHR1IqfT3OinW0VkjYisEpEF9ax7HxAa66c67S4VERUR710Wq6p+/UPttORbge5AGLAa6HdMm9uAZ12PJwBvOJ3bR/spDRgAvAJc5nRmH++rsUAb1+Of2e9Ug/3Uts7ji4BPnc7ti/3kahcDzAMWARneyhcII4nhQLaq5qhqBTATGH9Mm/HAy67HbwFnioh4MaMvaLSfVHW7qn4H1DgR0Ie401dfq2qp6+kialdiDDTu9NPhOk+jgEC8ksad7yiAPwKPAGXeDBcIRSIJyK/zvMD1Wr1tVLUKOATEeyWd73Cnn0ytpvbVjcAnHk3km9zqJxH5uYhsBf4G3OmlbL6k0X4SkSFAiqp+5M1gEBhFwhjHiMgkIAN41OksvkpVn1LVdOBe4H6n8/gaEQkC/gnc48TnB0KR2AGk1Hme7Hqt3jYiEgLEAoVeSec73OknU8utvhKRs4DfAheparmXsvmSpv5OzQQu9mgi39RYP8UA/YEsEdkOjARme+vkdSAUiaVATxHpJiJh1J6Ynn1Mm9nAda7HlwFfqetMUQBxp59MrUb7SkQGA89RWyD2OpDRF7jTTz3rPD0f2OLFfL7iuP2kqodUtYOqpqlqGrXnuC5S1WXeCOf3RcJ1juF2YA6wAZilqutE5CERucjV7AUgXkSygbuBBi9B81fu9JOIDBORAmrXIX9ORNY5l9g5bv5OPQpEA2+6Lu8MuILrZj/dLiLrRGQVtX971zWwO7/lZj85xqblMMYY0yC/H0kYY4xpPisSxhhjGmRFwhhjTIOsSBhjjGmQFQljjDENsiJhjDGmQVYkjDHGNMiKhDEtzHUT1P/Vef5nEbnLyUzGNJfdTGdMC3MtWvWOqg5xTc62BRiuqoE2H5jxAyFOBzDG36jqdhEpdM3flAistAJhWisrEsZ4xn+AyUAnYJqzUYxpPjvcZIwHuGbzXAOEAj1VtdrhSMY0i40kjPEAVa0Qka+Bg1YgTGtmRcIYD3CdsB5J7bTqxrRadgmsMS1MRPoB2cCXqhqIi+gYP2LnJIwxxjTIRhLGGGMaZEXCGGNMg6xIGGOMaZAVCWOMMQ2yImGMMaZB/x9I2qo7aoIVQwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "try : \n",
    "    # a sketch of the inflow profile:\n",
    "    import matplotlib\n",
    "    import numpy as np\n",
    "    import matplotlib.pyplot as plt\n",
    "    %matplotlib inline  \n",
    "\n",
    "    s = np.arange(0.0, 0.42, 0.01)\n",
    "    bvs = 6*s*(0.41-s)/(0.41)**2 \n",
    "    plt.plot(s, bvs)\n",
    "\n",
    "    plt.xlabel('y')\n",
    "    plt.ylabel('$u_x(0,y)$')\n",
    "    plt.title('inflow profile')\n",
    "    plt.grid(True)\n",
    "    plt.show()\n",
    "except ImportError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To solve the Stokes (and later the Navier-Stokes) problem, we introduce the bilinear form:\n",
    "\n",
    "$$\n",
    "  a((u,p),(v,q)) := \\int_{\\Omega} \\nu \\nabla u : \\nabla v - \\operatorname{div}(v) p - \\operatorname{div}(u) q\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "u,p = X.TrialFunction()\n",
    "v,q = X.TestFunction()\n",
    "\n",
    "a = BilinearForm(X)\n",
    "stokes = nu*InnerProduct(grad(u),grad(v))-div(u)*q-div(v)*p\n",
    "a += SymbolicBFI(stokes)\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(X)   \n",
    "f.Assemble()\n",
    "\n",
    "inv_stokes = a.mat.Inverse(X.FreeDofs())\n",
    "\n",
    "res = f.vec.CreateVector()\n",
    "res.data = f.vec - a.mat*gfu.vec\n",
    "gfu.vec.data += inv_stokes * res\n",
    "\n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## IMEX time discretization\n",
    "For the time integration we consider an semi-implicit Euler method (IMEX) where the convection is treated only explicitly and the remaing part implicitly:\n",
    "\n",
    "Find $(\\mathbf{u}^{n+1},p^{n+1}) \\in X = \\mathbf{V} \\times Q$, s.t. for all $(\\mathbf{v},q) \\in X = \\mathbf{V} \\times Q$\n",
    "\\begin{equation}\n",
    "   \\underbrace{m(\\mathbf{u}^{n+1},\\mathbf{v}) + \\Delta t ~\\cdot~a((\\mathbf{u}^{n+1},p^{n+1}),(\\mathbf{v},q))}_{ \\to M^\\ast} \n",
    "   ~=~ m(\\mathbf{u}^{n},\\mathbf{v}) - \\Delta t ~\\cdot~c(\\mathbf{u}^{n}; \\mathbf{u}^{n},\\mathbf{v}) \n",
    "\\end{equation}\n",
    "\n",
    "with \n",
    "\\begin{equation}\n",
    "   m(\\mathbf{u},\\mathbf{v}) = \\int \\mathbf{u} \\cdot \\mathbf{v}\n",
    "\\end{equation}\n",
    "and \n",
    "\\begin{equation}\n",
    "   c(\\mathbf{w},\\mathbf{u},\\mathbf{v}) = \\int \\mathbf{w} \\cdot \\nabla \\mathbf{u} \\cdot \\mathbf{v}\n",
    "\\end{equation}\n",
    "\n",
    "We prefer the incremental form (as it homogenizes the boundary conditions):\n",
    "\\begin{equation}\n",
    "   m(\\Delta \\mathbf{u}^{n+1},\\mathbf{v}) + \\Delta t ~\\cdot~a((\\Delta \\mathbf{u}^{n+1},p^{n+1}),(\\mathbf{v},q))\n",
    "   ~=~ \\Delta t (-a((\\mathbf{u}^{n},p^n),(\\mathbf{v},q)) -c(\\mathbf{u}^{n}; \\mathbf{u}^{n},\\mathbf{v}))\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dt = 0.001\n",
    "# matrix for implicit Euler \n",
    "mstar = BilinearForm(X)\n",
    "mstar += SymbolicBFI(InnerProduct(u,v) + dt*stokes)\n",
    "mstar.Assemble()\n",
    "inv = mstar.mat.Inverse(X.FreeDofs())\n",
    "\n",
    "velocity = gfu.components[0]\n",
    "absvelocity = Norm(velocity)\n",
    "\n",
    "conv = LinearForm(X)\n",
    "conv += SymbolicLFI(InnerProduct(grad(gfu.components[0])*velocity,v))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "t = 0\n",
    "tend = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t= 0.999000000000000876"
     ]
    }
   ],
   "source": [
    "# implicit Euler/explicit Euler splitting method:\n",
    "tend += 1\n",
    "while t < tend-0.5*dt:\n",
    "    print (\"\\rt=\", t, end=\"\")\n",
    "\n",
    "    conv.Assemble()\n",
    "    res.data = a.mat * gfu.vec + conv.vec\n",
    "    gfu.vec.data -= dt * inv * res    \n",
    "\n",
    "    t = t + dt\n",
    "    Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary 1: Computing stresses\n",
    "\n",
    "In the previously used benchmark the stresses on the obstacle are evaluated, the so-called drag and lift coefficients. \n",
    "\n",
    "The force acting on the obstacle is\n",
    "\n",
    "$$\n",
    "  F_{\\circ} = (F_D,F_L) = \\int_{\\Gamma_{cyl}} \\sigma_n = \\int_{\\Gamma_{cyl}} (-\\nu \\frac{\\partial u}{\\partial n} + p I) \\cdot n\n",
    "$$\n",
    "\n",
    "The drag/lift coefficients are \n",
    "\n",
    "$$\n",
    "  c_D = \\frac{2 }{\\bar{u} L} F_D, \\quad c_L = \\frac{2 }{\\bar{u} L} F_L.\n",
    "$$\n",
    "\n",
    "where $\\bar{u} = 1$ is the mean inflow velocity and $L = 0.41$ is the channel width."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We use the residual of our discretization to compute the forces. On the boundary degrees of freedoms of the disk we overwrote the equation by prescribing the (homogeneous) boundary conditions. The equations related to these dofs describe the force (im)balance at this boundary. \n",
    "\n",
    "Testing the residual (functional) with the characteristic function on that boundary (in the x- or y-direction) we obtain the integrated stresses (in the x- or y-direction):\n",
    "\n",
    "We define the functions which are characteristic functions (in terms of boundary evaluations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "drag_x_test = GridFunction(X)\n",
    "drag_x_test.components[0].Set(CoefficientFunction((-20.0,0)), definedon=mesh.Boundaries(\"cyl\"))\n",
    "drag_y_test = GridFunction(X)\n",
    "drag_y_test.components[0].Set(CoefficientFunction((0,-20.0)), definedon=mesh.Boundaries(\"cyl\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We will collect drag and lift forces over time and therefore create empty arrays\n",
    "\n",
    "and reset initial data (and time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "time_vals = []\n",
    "drag_x_vals = []\n",
    "drag_y_vals = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# restoring initial data\n",
    "res.data = f.vec - a.mat*gfu.vec\n",
    "gfu.vec.data += inv_stokes * res\n",
    "\n",
    "t = 0\n",
    "tend = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "With the same discretization as before."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Remarks:\n",
    " * you can call the following block several times to continue a simulation and collect more data\n",
    " * note that you can reset the array without rerunning the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t= 0.999000000000000876"
     ]
    }
   ],
   "source": [
    "# implicit Euler/explicit Euler splitting method:\n",
    "tend += 1\n",
    "while t < tend-0.5*dt:\n",
    "    print (\"\\rt=\", t, end=\"\")\n",
    "\n",
    "    conv.Assemble()\n",
    "    res.data = a.mat * gfu.vec + conv.vec\n",
    "    gfu.vec.data -= dt * inv * res    \n",
    "\n",
    "    t = t + dt\n",
    "    Redraw()\n",
    "   \n",
    "    time_vals.append( t )\n",
    "    drag_x_vals.append(InnerProduct(res, drag_x_test.vec) )\n",
    "    drag_y_vals.append(InnerProduct(res, drag_y_test.vec) )\n",
    "    #print(drag)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    " Below you can plot the gathered drag and lift values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xmc3XV97/HX5ywzk9kSkkmGJASCLBZBJSRsesVEbC9NFbyWClw3vNa0rrjU1to+lHrt417b6sMqKuJ1qUsJaqtNEeqaGLSATSBAAkgCJDBk32bmZJazfe4fv98MJ7OeOTm/czLzez8fj/M4v/WczzfLeZ/vb/kec3dEREQAEvUuQERETh4KBRERGaZQEBGRYQoFEREZplAQEZFhCgURERmmUBApk5l9w8w+We86RKKkUBARkWEKBZETZGapetcgUi0KBZFxmNkyM3vAzHrN7A6gKVy+0sy6zOwvzGwv8HUzO8XM7jSzA2Z2JJw+reS1zjSzjeFr/czMvmBm365X20TGo1AQGYOZNQA/BL4FzAW+B/xhySanhsvPANYQ/F/6ejh/OtAP3FKy/T8DvwHmATcDb460ASIVMo19JDKamV0BrAUWe/ifxMz+E/gF8DPgJ0C7uw+Ms/+FwHp3P8XMTgeeCrfvC9d/G8Dd3xR5Y0SmQD0FkbEtAp7z47817SqZPlAaCGbWbGZfNrNdZtYDbATmmFkyfK3DQ4EQejbK4kUqpVAQGdseYLGZWcmy00umR3axPwS8ELjU3duBK8LlFr7WXDNrLtl+SZXrFakKhYLI2O4F8sD7zCxtZq8HLplg+zaC8whHzWwu8PGhFe6+C9gE3GxmDWZ2OfDa6EoXqZxCQWQM7p4FXg/cCBwGrgP+dYJdPgvMAg4C9wH/MWL9G4HLgUPAJ4E7gMGqFi1SBTrRLFIH4SWuj7v7xyfdWKSG1FMQqQEzu9jMzjKzhJldBVxDcMmryElFd2KK1MapBIef5gFdwDvd/cH6liQymg4fiYjIMB0+EhGRYdPu8FFHR4cvXbq0on2PHTtGS0tLdQs6yanN8aA2x8OJtHnz5s0H3X3+ZNtNu1BYunQpmzZtqmjfDRs2sHLlyuoWdJJTm+NBbY6HE2mzme2afCsdPhIRkRIKBRERGaZQEBGRYQoFEREZplAQEZFhCgURERmmUBARkWHT7j6F6ay7L8f9Tx9i56FjDOSKzJ6V5uwFrbz4tNm0N6XrXZ6IiEKhFvZ09/PpnzzBui27yRaKo9ank8bLz+7gfyxbzOoXLySdVAdOROpDoRCx/9i6hw9/72FyxSLXX7KE17xkEb+zsI1Z6SRH+rL8dm8v92w/yI8e3sNNa7fwqbsf512rzub6i5eQUjiISI0pFCJ0+2+e4aM/eIQLl8zhs9ddyBnzjh+zZEFbEwvamnjFOfP5yFW/w/rf7udLG57kr3+4lW/du4ubrz6fy8+aV6fqRSSO9FU0Ij97dB9/9YNHWHnufG5/x2WjAmGkRMK48rxOvvenl3Prmy6iL5fnhq/cx83rttGfLdSoahGJO4VCBJ493Mf779jCBYtn84U3XkRTOln2vmbGVRcs5CfvfyU3vmwp3/jPnaz+3D08+MyRCCsWEQkoFKqsWHQ+/P2HAPjiGy+iuaGyI3SzGpLcfPX5/PM7LiWbL3Ltrffy+Z9vp1DUjyKJSHQUClX2wy3Pcd9Th/nrPziP005pPuHXe9lZHdz9/lfwmpcs5NM/fYLrb7uXZw/3VaFSEZHRFApV1J8t8Pc//i0vOW02b1ixpGqv296U5h+vX8Znr7uQx/b0svof7+HftjxXtdcXERmiUKii79y/iz3dA3x09XkkElb113/dssXcfdMrOPfUNm5au4UP3LGFnoFc1d9HROJLoVAl2XyRr/7qaS49cy6XvSC6y0iXzG3mjjWX8YFXn8u6h3az+h/vYdPOw5G9n4jEi0KhSv79od3s6R7gT1eeFfl7pZIJbnr1OXz3Ty7HDN7w5Xv5zE+fID/G3dIiIlOhUKiSb9+/i7MXtLLy3El/F7tqlp9xCne97xW8btliPvfz7bz+S/+pXoOInBCFQhVs39fLg88c5fqLl2BW/XMJE2lrSvOZN1zI529Yxt7uAa699V7e+e3NPL63p6Z1iMjMoGEuquC7m54llTBet2xx3Wp47UsXceV5C/jKxqf58sYnuXvrXq44dz6XtBe4ouiRnPgWkZlHoXCCcoUiP3jwOa48bwEdrY11raW5IcVNrz6Ht77sDL5z/zN8/dc72fjEILfvWM81Fy7i6gsX8cLOtpr3ZkRk+ogsFMysCdgINIbv8313//iIbRqBbwLLgUPAde6+M6qaonD/U4c5mMny+otOq3cpw+Y0N/DuVWfz9v92Jp/93noeH2zlyxuf4osbnmTR7CZe+cL5XHrmPF66ZA5L5zUrJERkWJQ9hUHgVe6eMbM08Cszu9vd7yvZ5u3AEXc/28yuBz4FXBdhTVX34217mZVO8soanmAuV1M6yWWLUnxk5SUc6B3kF4/vY/3jB7jzoT3c/ptnAWhvSnH+otmctaCFMztaeUFHC0s7Wlg4u2lKYzaJyPHcnZ6BPD39OfJFp1B03J2CB9PFIphBKmmkEkYykSCVMFJJI5kwUolE+GzhNrU5BRxZKLi7A5lwNh0+Rg7ccw1wczj9feAWM7Nw35Nesej89NF9XHFux0n/ATq/rZHrLj6d6y4+nXyhyPb9GR7uOsqWZ7t5bE8P67bspmcgf9w+c5rTnNreRGd7E53tjcxva2RuSyNzW9Kc0tzA3JaG4efmhqR6HDIjuTvHsgWO9mU52pfjaF+OI31ZjvbnOHoseD7Sl6W7dHlfju7+XNXHKlt9ZpqVK6v6kqNEek7BzJLAZuBs4Avufv+ITRYDzwK4e97MuoF5wMEo66qWh5/rZm/PAH9+/gvrXcqUpJIJzlvYznkL27nu4mCZu3OkL8dTBzLsPNTH3u5+9vYMsLd7kH09Azy6p4dDmUHG+zfemEowt6WBOc0NtDWlaG9K0d6Upq0pRduI5/ZZ6eFtWhvTzEonaWpI0JBMKFgECL5w9ecK9OcK5AtOrlAMv20XyRWcfMHZcbRAy87D4bfukm/g7hSKlEwHz0PLR257bLBA70COnoEcPf15egZy9Ibf8I/05ejuz5IrjP/h3tKQZE5zA3Oagy9LC+fM4pRwevasNLNnpUknE5hBMmEkzTALegMe1jHUk8gVisfN54tOPmx7vuCkjj4T+Z+91eJLuZnNAX4AvNfdt5Ys3wpc5e5d4fyTwKXufnDE/muANQCdnZ3L165dW1EdmUyG1tbWyhoxhu8/keWup3N8blUzrQ0n54dZNdtcdKc/D71ZpzfrZHLhc9bpzRE8Z53+vNOXh/58OJ0b3UUciwENyfCRMBqTwU+VNobz6SSkEpAySCWMdCKcT1j4DOmEUcgN0tLUODyfKtnu+X3GeA17fl1imoVTtf9tV6JQDP59DP29D00P5KEv7wyULBssPP88WAi2GSwE6wcLwfJaSho0p6A5bTSnjFlpaE4ZLWmjNW20NEDr0HTaaG0wWtLBslQNr+w7kb/nVatWbXb3FZNtV5Orj9z9qJmtB64Ctpaseg5YAnSZWQqYTXDCeeT+twG3AaxYscJXVth/2rBhA5XuO5ZPP/IrVpyR5DW/d3nVXrPaqt3mSrg7fdlC8O1rIBd+K8vTO5AnM5BnIPxGOBA++nMF+rNFBvIFBrKF4W+Mg7kivYUi2XyRXPiczRcZLBTI5kvv5jYge0I1p5NGQzJBOhX0YBpS4SM54jkVHAc2MxIWhEnCDDOOWzYqY3zCWUZ+WfPj1o3ed/++ARZ0zp5037Hf+/kFpbuWThfcGcwXGcwVGMwXGcgFf+ZD08eyeQZyk99RnzBoaUzR2piipTFFS1OSeQ0pWhqTNDeEyxqSNIfPTekk6WRi+Lh7KpkIgj6Z4LFtj7DswpeSNCORCL55J8Jv4MGy0d/MRy5PJIyWhhRN6enRS63F/+corz6aD+TCQJgF/C7BieRS64C3AvcC1wK/mC7nE472Zdm6u5v3X3luvUs56ZlZ8J+9McWps5sieQ/3oKudzRdZ/8t7uPiyy4PQKAmP4enSZeH8UMgMli4bse3gGMsyg3nyheBj1X3oMEXQq3IfWhZ8ix75mTNqHptkfem641f29xU5kO8ee+PRs6P2P/61R9dkBo3pJE2pBO2z0sxva6QxlaApnaQhlaClIUlbU5rWxhStTSnaGoPDha1NQQAEhw9TzEpX79xTev9jvOKck+8Cj+kuyp7CQuCfwvMKCeC77n6nmX0C2OTu64CvAt8ysx3AYeD6COupqnufPIQ7vPxs/YbyycDMSCeNdDJBa4PR2R5N+JysToYeocwMUV599DCwbIzlHyuZHgD+KKoaovTrJw/S0pDkpUvm1LsUEZGq0dhHFfr1jkNc+oJ5pJP6IxSRmUOfaBXY093P0weP8bKzdOhIRGYWhUIFNu86AsAlZ86tcyUiItWlUKjAA7uO0pQObgATEZlJFAoVeOCZI7xk8RydTxCRGUefalM0kCuwbXc3y87QVUciMvMoFKZo2+5ucgXnotNPqXcpIiJVp1CYogd2HQVQKIjIjKRQmKItzx7ltFNmMb+tvr+yJiISBYXCFG3d3c2LF8+udxkiIpFQKExBz0COXYf6OH+RLkUVkZlJoTAFj+3uAeD8ReopiMjMpFCYgm3DoaCegojMTAqFKdi2u4eO1kYWxGxYZhGJD4XCFGzb3a1egojMaAqFMg3mC+zYn1EoiMiMplAo0/Z9GfJF10lmEZnRFApl2r6/F4AXntpa50pERKKjUCjT9n0ZUgnjjHkt9S5FRCQyCoUybd+f4cyOFg2XLSIzmj7hyrR9Xy/ndOrQkYjMbAqFMgzkCjxzuI+zF7TVuxQRkUgpFMrw1IFjFB3OWaCegojMbAqFMgxdeXRup3oKIjKzKRTKsGN/hmTCWNrRXO9SREQipVAowxP7ejljXjONqWS9SxERiZRCoQzb92d0PkFEYkGhMIlsvsiuQ32crVAQkRhQKEyi60gfhaLzgg6FgojMfJGFgpktMbP1ZvaomW0zs5vG2GalmXWb2Zbw8bGo6qnUzkPHAHSSWURiIRXha+eBD7n7A2bWBmw2s5+6+6MjtrvH3V8TYR0nZOfBPgCWaswjEYmByHoK7r7H3R8Ip3uBx4DFUb1fVHYeOkZbY4q5LQ31LkVEJHLm7tG/idlSYCNwgbv3lCxfCfwL0AXsBv7M3beNsf8aYA1AZ2fn8rVr11ZURyaTobV1aucG/mHTAJmsc/PLZlX0nvVWSZunO7U5HtTmqVm1atVmd18x6YbuHukDaAU2A68fY1070BpOrwa2T/Z6y5cv90qtX79+yvtc8Xe/8Hd/Z3PF71lvlbR5ulOb40Ftnhpgk5fxmR3p1UdmliboCXzH3f91jEDqcfdMOH0XkDazjihrmopcoUjXkX7O7ND5BBGJhyivPjLgq8Bj7v6ZcbY5NdwOM7skrOdQVDVNVdeRfgpF1w/riEhsRHn10cuBNwOPmNmWcNlHgdMB3P1W4FrgnWaWB/qB68Nuzklh58HgctQzdTmqiMREZKHg7r8CbJJtbgFuiaqGEzV0j4J6CiISF7qjeQI7DwaXo87T5agiEhMKhQnsPNTHGR3NhKc9RERmPIXCBHYdOqZDRyISKwqFcRSLznNH+1lyik4yi0h8KBTGsb93kFzBOe2U6Xkns4hIJRQK4+g6EgyEp1AQkThRKIyj60g/AEvm6vCRiMSHQmEczx4OegqL56inICLxoVAYR9eRfua3NdKUTta7FBGRmlEojKPraJ/OJ4hI7CgUxtF1pJ/TdDmqiMSMQmEMhaKz+2i/egoiEjsKhTHs6xnQPQoiEksKhTEMXY6qw0ciEjcKhTEM3bi2RD0FEYkZhcIYhnoKi3SPgojEjEJhDF1H+ligexREJIYUCmN47mi/egkiEksKhTHs6R5g0ZymepchIlJzCoUR3J09Rwc4tV09BRGJH4XCCN39OfpzBfUURCSWFAoj7OkeAGDhbPUURCR+FAoj7OkOLkc9dbZ6CiISPwqFEXYfDXoKOnwkInGkUBhhb/cAyYSxoE2hICLxo1AYYXd3PwvaGkkmrN6liIjUnEJhhL3dAyzU+QQRiSmFwgh7ugdYqLuZRSSmygoFM5s7xiM9yT5LzGy9mT1qZtvM7KYxtjEz+5yZ7TCzh83sokobUg3uwY/rLGxXT0FE4qncnsIDwAHgCWB7OL3TzB4ws+Xj7JMHPuTuLwIuA95tZi8asc3vA+eEjzXAl6ZYf1Ud7csxmC+qpyAisVVuKPwUWO3uHe4+j+DD/E7gXcAXx9rB3fe4+wPhdC/wGLB4xGbXAN/0wH3AHDNbWEE7qmJ3eI/CIp1TEJGYKjcULnP3Hw/NuPtPgMvDD/LGyXY2s6XAMuD+EasWA8+WzHcxOjhqZm94N7NuXBORuEqVud0eM/sLYG04fx2wz8ySQHGiHc2sFfgX4P3u3lNJkWa2huDwEp2dnWzYsKGSlyGTyUy47y+fyQGw89EH6X5qZpyDn6zNM5HaHA9qc0TcfdIH0AF8HngwfNwCzAcagLMn2C8N/Bj44DjrvwzcUDL/W2DhRLUsX77cK7V+/foJ13/q7sf8rL/8kecLxYrf42QzWZtnIrU5HtTmqQE2eRmf92X1FNz9IPDecVbvGGuhmRnwVeAxd//MOPuuA95jZmuBS4Fud99TTk1R2Ns9QGd7k25cE5HYKisUzGw+8OfA+cDwAXd3f9UEu70ceDPwiJltCZd9FDg93PdW4C5gNUGw9AFvm2L9VbW7u183rolIrJV7TuE7wB3Aa4A/Bd5KcFnquNz9V8CEX7nDLs27y6whcvt6BnnRovZ6lyEiUjflnk2d5+5fBXLu/kt3/1/ARL2EaWlfzwCdGghPRGKs3J5CLnzeY2Z/AOwG5kZTUn1kBvP0ZQt0tk96ha2IyIxVbih80sxmAx8iuAqpHfhAZFXVwb6e4B6FBQoFEYmxSUMhvBfhHHe/E+gGVkVeVR3s7xkE0OEjEYm1Sc8puHsBuKEGtdTV/t6hnoJCQUTiq9zDR782s1sIrkA6NrTQw7GNZgIdPhIRKT8ULgyf/yZ8NsCZQVcg7e8ZZFY6SVtjuX8kIiIzz4SfgGb2wXDyToIQKL3vwKMqqh729Q7S2d5IcCO2iEg8Tfa1uC18fiFwMfBvBMHwWuA3EdZVc/t7Bligk8wiEnMThoK7/w2AmW0ELvLgdxEws5uBH0VeXQ3t7x3kfN3NLCIxV+4dzZ1AtmQ+Gy6bMfappyAiUvaJ5m8CvzGzH4TzrwO+EUlFdaC7mUVEAuUOnf23ZnY38Ipw0dvc/cHoyqqtoctRO3WPgojEXNnXX4b3JMyY+xJKDd3NvKBNPQURibeZ8ZuTJ0h3M4uIBBQKlB4+Uk9BROJNocDzdzO36m5mEYk5hQK6m1lEZIhCAd2jICIyRKEAHOgd1OioIiIoFIDwt5l15ZGIiEJh6G5m3aMgIqJQ0N3MIiIlYh8KuptZROR5sQ+FA5kgFOYrFEREFAoHe4NQ6GhVKIiIKBQyg6QSxuxZ6XqXIiJSdwqFzCDzWhtIJHQ3s4hI7EPhQO+gzieIiIQiCwUz+5qZ7TezreOsX2lm3Wa2JXx8LKpaJnIwk9X5BBGRUJQ9hW8AV02yzT3ufmH4+ESEtYzrYGZQoSAiEoosFNx9I3A4qtevBnfnYEaHj0REhpi7R/fiZkuBO939gjHWrQT+BegCdgN/5u7bxnmdNcAagM7OzuVr166tqJ5MJkNra+vz81nnPb/o44bfaeC/L52ZVx+NbHMcqM3xoDZPzapVqza7+4pJN3T3yB7AUmDrOOvagdZwejWwvZzXXL58uVdq/fr1x81v39fjZ/zFnf7DB7sqfs2T3cg2x4HaHA9q89QAm7yMz9i6XX3k7j3ungmn7wLSZtZRyxoO9GYBmK9zCiIiQB0vSTWzUy38qTMzuySs5VAta9AQFyIix4vsR4nN7HZgJdBhZl3Ax4E0gLvfClwLvNPM8kA/cH3YxakZDXEhInK8yELB3W+YZP0twC1RvX85NMSFiMjxYn1H89A9ChriQkQkEOtQONA7SEdbQ73LEBE5acQ6FDTEhYjI8WIeChriQkSkVGxDwTXEhYjIKLENhe7+HLmCq6cgIlIitqFwMDN0j4JONIuIDIltKGiICxGR0eIbChriQkRklNiGgoa4EBEZLb6hoCEuRERGiXUoaIgLEZHjxTgUsszTlUciIseJbSgcOpZlns4niIgcJ76hkBlkXot6CiIipWIbCoePZRUKIiIjxDIU+rMF+rIF5uqcgojIcWIZCoeOBfcoqKcgInK8WIbC4WPBEBfzWnSiWUSkVCxD4VAmCAUdPhIROV48Q2G4p6BQEBEpFctQODx0TkH3KYiIHCeWoXAok6UhlaClIVnvUkRETirxDIXwHgUzjXskIlIqlqFw+JjGPRIRGUssQ+FQZpC5uhxVRGSUeIaChrgQERlTPEMho1AQERlLZKFgZl8zs/1mtnWc9WZmnzOzHWb2sJldFFUtpfqzBfpzGvdIRGQsUfYUvgFcNcH63wfOCR9rgC9FWMswjXskIjK+yELB3TcChyfY5Brgmx64D5hjZgujqmfI0BAXGvdIRGS0ep5TWAw8WzLfFS6L1NBgeDp8JCIyWqreBZTDzNYQHGKis7OTDRs2VPQ6mUyGLZsfAmDH1gfpeWrmn2fPZDIV/3lNV2pzPKjN0ahnKDwHLCmZPy1cNoq73wbcBrBixQpfuXJlRW+4YcMG5rctgUceZ/WVV9DaOC0y8YRs2LCBSv+8piu1OR7U5mjU86vyOuAt4VVIlwHd7r4n6jc9fEzjHomIjCeyr8pmdjuwEugwsy7g40AawN1vBe4CVgM7gD7gbVHVUupgJkuHxj0SERlTZKHg7jdMst6Bd0f1/uM5fGxQJ5lFRMYx88+0jnD4WFbjHomIjCN2oTB0+EhEREaLXSgEPQWFgojIWGIVCoN517hHIiITiFUo9GQdgA6dUxARGVOsQqE3F4SCDh+JiIwtVqHQMxiGgg4fiYiMKVah0KvDRyIiE4pXKOTUUxARmUisQqFnEI17JCIygViFQm/WNe6RiMgEYhcKOnQkIjK+2IWCfoZTRGR8sQqFnqwzT/coiIiMK1ah0Jt13bgmIjKB2ITCQK5AtginKBRERMYVq1AAmJXW5agiIuOJTSj0h6HQpFAQERlXbEJhIFcEoCkdmyaLiExZbD4hB9RTEBGZVOxCQecURETGF6NQCA4fNerwkYjIuGLzCTmQ1+EjEZHJxCYUBofOKaQUCiIi44lNKMxva2RFZ1J3NIuITCA2obD8jLm8Z1kTp85uqncpIiInrdiEgoiITE6hICIiwxQKIiIyLNJQMLOrzOy3ZrbDzD4yxvobzeyAmW0JH38cZT0iIjKxVFQvbGZJ4AvA7wJdwH+Z2Tp3f3TEpne4+3uiqkNERMoXZU/hEmCHuz/l7llgLXBNhO8nIiInyNw9mhc2uxa4yt3/OJx/M3Bpaa/AzG4E/g9wAHgC+IC7PzvGa60B1gB0dnYuX7t2bUU1ZTIZWltbK9p3ulKb40FtjocTafOqVas2u/uKybaL7PBRmf4duN3dB83sT4B/Al41ciN3vw24DWDFihW+cuXKit5sw4YNVLrvdKU2x4PaHA+1aHOUofAcsKRk/rRw2TB3P1Qy+/+Av5vsRTdv3nzQzHZVWFMHcLDCfacrtTke1OZ4OJE2n1HORlGGwn8B55jZmQRhcD3wP0s3MLOF7r4nnL0aeGyyF3X3+ZUWZGabyuk+zSRqczyozfFQizZHFgrunjez9wA/BpLA19x9m5l9Atjk7uuA95nZ1UAeOAzcGFU9IiIyuUjPKbj7XcBdI5Z9rGT6L4G/jLIGEREpX9zuaL6t3gXUgdocD2pzPETe5sguSRURkeknbj0FERGZgEJBRESGzchQKGMgvkYzuyNcf7+ZLa19ldVVRps/aGaPmtnDZvZzMyvrmuWT2WRtLtnuD83MzWzaX75YTpvN7A3h3/U2M/vnWtdYbWX82z7dzNab2YPhv+/V9aizWszsa2a238y2jrPezOxz4Z/Hw2Z2UVULcPcZ9SC4/PVJ4AVAA/AQ8KIR27wLuDWcvp5gUL661x5xm1cBzeH0O+PQ5nC7NmAjcB+wot511+Dv+RzgQeCUcH5BveuuQZtvA94ZTr8I2Fnvuk+wzVcAFwFbx1m/GrgbMOAy4P5qvv9M7CmUMxDfNQRDagB8H7jSzKyGNVbbpG129/Xu3hfO3kdwh/l0Vu6Ai/8b+BQwUMviIlJOm98BfMHdjwC4+/4a11ht5bTZgfZwejawu4b1VZ27byS4b2s81wDf9MB9wBwzW1it95+JobAYKB1UrytcNuY27p4HuoF5NakuGuW0udTbCb5pTGeTtjnsVi9x9x/VsrAIlfP3fC5wrpn92szuM7OralZdNMpp883Am8ysi+C+qPfWprS6mer/9ymp94B4UmNm9iZgBfDKetcSJTNLAJ8hfnfJpwgOIa0k6A1uNLMXu/vRulYVrRuAb7j7p83scuBbZnaBuxfrXdh0NBN7CpMOxFe6jZmlCLqch5i+ymkzZvZq4K+Aq919sEa1RWWyNrcBFwAbzGwnwbHXddP8ZHM5f89dwDp3z7n70wRD0p9To/qiUE6b3w58F8Dd7wWaCAaOm6nK+v9eqZkYCsMD8ZlZA8GJ5HUjtlkHvDWcvhb4hYdncKapSdtsZsuALxMEwnQ/zgyTtNndu929w92XuvtSgvMoV7v7pvqUWxXl/Nv+IUEvATPrIDic9FQti6yyctr8DHAlgJmdRxAKB2paZW2tA94SXoV0GdDtzw8sesJm3OEjL28gvq8SdDF3EJzQub5+FZ+4Mtv890Ar8L3wnPoz7n513Yo+QWW2eUYps80/Bn7PzB4FCsCH/fgh6qeVMtv8IeArZvYBgpPON07nL3lmdjtBsHeE50k+DqQB3P1WgvMmq4EdQB/wtqq+/zT+sxMRkSqbiYePRESkQgoFEREZplAQEZFhCgVOmST+AAABRUlEQVQRERmmUBARkWEKBZFJmNkcM3tXOL3IzL5f75pEoqJLUkUmEQ6tfqe7X1DnUkQiN+NuXhOJwP8FzjKzLcB24Dx3v8DMbgReB7QQDCXxDwTDO78ZGARWu/thMzsL+AIwn+Bmo3e4++O1b4bI5HT4SGRyHwGedPcLgQ+PWHcB8HrgYuBvgT53XwbcC7wl3OY24L3uvhz4M+CLNalapALqKYicmPXu3gv0mlk38O/h8keAl5hZK/Aynh9eBKCx9mWKlEehIHJiSkebLZbMFwn+fyWAo2EvQ+Skp8NHIpPrJRiKe8rcvQd42sz+CIZ/X/el1SxOpJoUCiKTCEcZ/XX4Q+p/X8FLvBF4u5k9BGxj7J8NFTkp6JJUEREZpp6CiIgMUyiIiMgwhYKIiAxTKIiIyDCFgoiIDFMoiIjIMIWCiIgM+//nTAtIcyuohwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEWCAYAAACaBstRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8W+d18PHfA3DvTXGJlMQhUVuiliXZsi3POlIcN6vZTeo2q23SpEmat3nbJG1WR5o36XDSNquNE2c4tmNbtmXRsmxLlixSk6IoUVyiuPcmgOf9A4BM06QIksC9F8D5fj78CCQv7z1XIHHwrPMorTVCCCGEL2xmByCEECJ4SNIQQgjhM0kaQgghfCZJQwghhM8kaQghhPCZJA0hhBA+k6QhhJ8ppRqUUnuVUn+llPrBlK/fr5RqVkoNKaU2mhmjEAslSUOIANFa/73W+iNTvvQPwCe01gla6ypvcjErPiEWQpKGEMYpBM6ZHYQQiyFJQ4gAUUr9jVLqp0qpaKXUEGAHTimlLiulfgIsBR73dFf9pbnRCuGbCLMDECLUaa3HgQSllAbWa60vASildgMf0Vo/Z2qAQsyDtDSEEEL4TJKGEEIIn0nSEMI8UmJaBB1JGkKYpx1YbnYQQsyHJA0hzPM14P8opfqUUp8xOxghfKFkEyYhhBC+kpaGEEIIn0nSEEII4TNJGkIIIXwmSUMIIYTPQq6MSEZGhi4qKlrQzw4PDxMfH+/fgCxO7jk8yD2Hh8Xc82uvvdaltc6c67iQSxpFRUWcOHFiQT9bWVnJnj17/BuQxck9hwe55/CwmHtWSjX6cpx0TwkhhPCZJA0hhBA+k6QhhBDCZ5I0hBBC+MzUpKGUulspVauUuqSU+vwNjntAKaWVUhVGxieEEOKNTEsaSik78D3gHqAceLdSqnyG4xKBPwOOGRuhEEKI6cxsaWwFLmmt67XWE8DDwP4ZjvsK8A1gzMjghBBCvJmZSSMPaJ7yeYvna9cppTYBBVrr3xkZmBDCWD3DE/z8eBPffb6OA+famHS6zA5JzMKyi/uUUjbgn4AP+nDsg8CDANnZ2VRWVi7omkNDQwv+2WAl9xwerHrPWmueb3bw89oJJpyvfz07TvHR9dEUJdsXfG6r3nMgGXLPWmtTPoAdwIEpn38B+MKUz5OBLqDB8zEGtAIVNzrv5s2b9UIdOnRowT8brOSew4NV7/lbT1/QhZ97Qr//P4/p8639emTcoZ8516Zv+tpBXfZ/ntTHr3Qv+NxWvedAWsw9Aye0D6/dZnZPHQdKlFLLlFJRwLuAx7zf1Fr3a60ztNZFWusi4CiwT2u9sBohQghL+U1VC989dIl3bSngvz+4hVU5ScRG2bmjPJtHP76T3ORYPvLjE7QPyHCmlZiWNLTWDuATwAGgBviF1vqcUurLSql9ZsUlhAi8hq5hvvDrM2xblsZX37oGm0294fuZidF8/wMVjE06+cwjp7y9D8ICTF2nobV+UmtdqrVeobX+O8/XvqS1fmyGY/dIK0OI4Ke15kuPnSPCZuM7795IhH3ml6EVmQn81b2reLGui9+duWZwlGI2siJcCGGoZ8+3c/hiJ39xZynZSTE3PPY92wpZuSSRbz5dy7jDecNjhTEkaQghDKO15tvP1VGUHsf7thfOebzdpvjCvato6hnhkRMtBkQo5iJJQwhhmOdqOjh/bYBP3FYya7fUdDeXZLC+IIUfvFiP0yVjG2aTpCGEMMx3D11iaVocb92Q6/PPKKV4cPdyGrpHePZ8WwCjE76QpCGEMMSp5j5ONffxkd3LfG5leN21Opu8lFh+erQpQNEJX0nSEEIY4idHG4mLsnP/xry5D54mwm7j7RX5vHS5i5bekQBEJ3wlSUMIEXB9IxM8fqqV+zfmkRgTuaBzPLApH63hV69d9XN0Yj4kaczBJQNvQizao1VXGXe4eM+2uWdMzaYgLY6dxen88mSz/F2aSJLGDLTW/M+xRvZ86xDFX3ySvf/0Aj97tUl+UYVYoN+eamVVThLluUmLOs/vb86nuWeUquY+P0Um5kuSxjRaa77w6zN88TdnyUiI5sGbV5AQHcEXfn2GT/6sigmHlGwWYj6aukeoaupj/zxmTM3m9lXZRNoVT8kKcdNYtjS6Wb7/Yj0PH2/mY3tW8Nm7ylBKobXm+y/W8/dPXiDSrvjnd25AKTX3yYQQPHbKPQbxlvWLTxpJMZHsKs7gqbNtfPH3VsnfoQmkpTFF27CLbx2o5a7V2dcTBnjmid+8gs/cWcqj1a08fLx5jjMJIcDdcv9tdStbilLJS4n1yznvWZPD1b5Rzl4d8Mv5xPxI0pji13UTxETY+cpb18z4DuZje4rZVZzBlx8/T1O3TPsTYi617YPUdQyxzw+tDK87yrOx2xRPnZUuKjNI0piia1SzYWkKWYkzF1Gz2RTfevs6AP7+yRojQxMiKD17rh2Au9Ys8ds5U+OjqChM5VBtp9/OKXwnSWMKDdjm6CPNSY7l47eu4OlzbRyt7zYmMCGC1HM17WwomP2N2ELtKcui5tqAbNBkAkkaU2nwZVztI7uXk5UYzXcO1gU+JiGCVPvAGKda+rmjPNvv576lNBOAwxeltWE0SRpTaMCXuRgxkXb+aPdyXr7czcmm3kCHJURQOljTARCQpLEqJ5HMxGhekKRhOEkaU2jweQrfH2xbSkpcJP966FJggxIiSD1X087StDhKshL8fm6lFLeUZvJiXZeUSzeYJI1pfJ31HR8dwft3FHHwQofMpBJimtEJJ0cudXH7qqyAraW4pTST/tFJTrXI6nAjSdKYQvs4puH17q0FKOBnx6VcsxBTHbvSzYTDxa1lWQG7xu6SDJSCFy92Bewa4s1kRfgU7kau71kjJzmW21Zm88iJZj61t5SoCMnBZmrtG+Xh480cre+mb2SCjIRoti1L511bC+bci1r414t1XURF2Ni6LC1g10iJi2LVkiSOXekGSgJ2HfFG8io3hdYa2zxb0u/ZtpSuoQmePd8emKDEnCadLv7xmVr2fKuS7z5fx7jDxbKMeIbGHXz74EV2f/MQ/1Z5Wfq+DXSkroutRWnERNoDep3ty9N5rbGXcYczoNcRr5OWxjTz7X69uTST7KRoflN1ld9blxOYoMSshsYdfORHxzla38PbNubx6TtLyU+Nu/79xu5h/v7JGr7x9AWqm3v5l3dtNDHa8NA+MEZt+yD3b5r/ZkvztW15Gv/10hVOt/SzpShwrRrxOmlpTOGecju/rGG3Kd6yLpcXLnbQNzIRmMDEjAbHJnn/fx7jeEMv//SO9fzTOze8IWEAFKbH8x/vq+Cv7yvnmfPtfPSnr+GQFkdAHalzjzHsLskI+LW2ehLFMVloaxhJGlO4p9zO/+f2b8hj0ql56qxsem8Ul0vz5w9Xc7qln+/9wSbetin/hsd/eNcyvvrWNRyq7eSnNZLcA+nFuk7S493jDYGWGh/FyiWJHK3vCfi1hJskjanmOXvKa01eEssz4vlttWxDaZR/fu4iBy908KW3lHO3j3WN3rOtkD++ZTmVzQ4eP9Ua4AjDk8ulOXKpm10lGdjmO0C4QN5xjUmn7HVjBEkaUyykewrcC43uW5/LsSs9dA2N+z8w8QYnGnr47qFLvH1zPu/bPr/tQz9zZxkrkm381a/PcLVvNEARhq+LHYN0DY2zqzjwXVNe25alMTrp5HRLv2HXDGeSNKbQMJ8Zt29wZ3k2WsPBGplFFUijE04+88gp8lJi+b/7Vs974Vik3cafrI9m0uXiy4+fC1CU4evoZffYwvbl6YZds8IzrnGyUUr6GEGSxhRaLzhnsDo3ibyUWJl6G2DfO3SJhu4RvvnAOhKiFzb5LzPOxp/eXsKBc+08f0GeL386dqWHvJRYCtLi5j7YTzIToylIi5U6cAaRpDHNXKXRZ6OU4o7ybF6s62JkwuHnqARAS+8I33+xnv0bcrlpkd0fH9m1nOKsBL76RA0O6Qv3C601r17pYVsAF/TNZmNBKiebetFaZsYFmiSNKRY6e8rrzvJsxh0uDktZg4D4xtO1KAWfu3vlos8VFWHjL+8qo75rmF+dbPFDdOJSxxDdwxNsW2580ti0NIX2gXGu9cv+GoEmSWOKxXRPAWxZlkZSTATPnJept/52qrmPx0+18uDu5eT6aa/pO8qz2VCQwr88V8fYpKwoXqyjV9zTXrctM248w2tTYSqAdFEZQJLGNIupyBlpt3HryiwOX+zEJQvI/Oo7B+tIiYvkwVtW+O2cSik+e1cZrf1j0trwg1ev9JCdFE1hunHjGV6rcpKIjrBxslEq3gaaJI0pfN2E6UZuLsmka2iC89cG/BGSAM5e7efghQ4+smvZgge/Z3PTinTW5yfzgxevSG2qRdBac6y+m23L0gNWCv1GIu021uUnS0vDAJI0ptB+yBre0gmH62RHMX/5l4N1JMVE8P6bivx+bqUUD968gitdwzLzbREaukfoGBw3ZTzDa9PSVM63DkjxwgAzNWkope5WStUqpS4ppT4/w/c/rZQ6r5Q6rZQ6qJSa30quhcS0yKyRlRTDyiWJUuPfTy51DPHs+XY+uHMZSTGRAbnG3WuWsDQtjocOXw7I+cOBt/aTGeMZXhuXpjLhdHH2qrTyA8m0pKGUsgPfA+4ByoF3K6XKpx1WBVRordcBvwS+GciYNMy7NPpMbinN5ERjD8PjMvV2sX748hWi7DbevyNw7xfsNsWHdy3jZFMfp5qlT3whXm3oISMhihWZ8abFsGlpCgBV0kUVUGa2NLYCl7TW9VrrCeBhYP/UA7TWh7TW3r1UjwI3rkq3SPPduW82u0symXRqz+YwYqH6Ryb51WtX2b8hl4yE6IBe622b8oiLsvO/x2QXxoWobupjQ0GqKeMZXllJMWQnRXP2qpQTCSQz99PIA5qnfN4CbLvB8R8GnprpG0qpB4EHAbKzs6msrFxQQC7toq2tjcrKxb1TmXBqomzwv89XY2sL7IvdYg0NDS34/yvQnroyyeikk7XRXX6NcbZ7rshS/KaqmVuSu4mLNO/FLxAC+TwPT2rqu0bYmDph+u9SboyDY3XXqKzst/TvdqAYcc9BsQmTUuq9QAVwy0zf11o/BDwEUFFRoffs2bOw6xx6ktycHPbsWbfASF93U9Or1HePsNBYjFJZWWnJGB1OF188Wsm2ZWm8f98Ov557tntOK+5j33dfojthGffuKPLrNc0WyOf5hYudwKu87ZZN7DSwUOFMTjnq+PbBi1Ts2MWJV45Y8nc7kIz4ezaze+oqUDDl83zP195AKbUX+CKwT2sd0BKyi10RPtWu4gzqu4ZpH5AVqgvxwsVOrvaN8qGdRYZdc11+CmvzkvmfY01SjmIeqpv6UArW5SebHQrr8pPRGs5JF1XAmJk0jgMlSqllSqko4F3AY1MPUEptBP4Dd8LoCHRAejFlbqfxziI5KjuKLcjPjzeTkRDN7auyDb3uO7YUcKFtkJprg4ZeN5hVN/dSkpVAYoBmt83Hmjx34jojSSNgTEsaWmsH8AngAFAD/EJrfU4p9WWl1D7PYd8CEoBHlFLVSqnHZjmdv6LyW0ujPDeJxOgI2VFsAToGx3j+QgcPbMoj0m7sr+h9a3OItCt+UyUrxH2htaa6uY8NBSlmhwK4K94uSYqRpBFApo5paK2fBJ6c9rUvTXm819B48Fc7wz2Nc+uyNJlBtQC/PnkVh0vz9oqCuQ/2s9T4KPaUZfHb6lY+f88q7AbtPhesGrtH6B2ZZENBqtmhXLc2P5kzV/u537cNHcU8yYrwKbReeGn0mWxfnk595zAdMq7hM601vzjeTEVhKsVZCabEcP/GPDoGx3n5sizQnEu1Z12LVVoaAGvzkqnvHGbUIeNSgSBJYwp/DoTD67uXeat/irmdaOylvmuYd2wxvpXhddvKLBJjIvjNSdnzfS7VzX3ERtopzTYnwc9krWdAvnFA9kkJBEka0/izM+L1cQ3povLVIyeaiY+y83trc0yLISbSff2nz7UxOiF1jG6kqrmPtfnJRBg89nQjaz2D4Q39kjQCwTrPtAW4V4T7L23YbYoty9Ikafho3OHkqbNt3LVmCfF+rmY7X29Zn8vIhJMXLgZ80l7QGnc4qWkdYONS63RNAWQkRJObHEPDgCT8QJCkMUUgekC3LUujvnOYzsGALjEJCZW1nQyOOdi/Ic/sUNi2LI3UuEieOisbas3mXOsAE04XGy00nuG1Ji9ZWhoBIkljGn+XzqkocpeKljr/c3usupX0+Ch2rjCvUqpXhN3GneVLOFjTIaW2Z1Hd5B0Et87MKa/y3CTaRzQjE1I01N8kaUzh3u7Vv1ljTV4SUXYbJxsladzI4Ngkz9W083vrcizTP37P2iUMjTs4UiezqGZS3dzHkqQYliTHmB3Km6zKSUIDF9pkkaa/WeOv00L83dKIjrCzNj+ZE5I0bujZ8+2MO1zs35BrdijX3bQig8SYCJ48I11UM7HSor7pynOSAKiRHTT9TpLGFC78s5/GdJsLUznT0i/dHDfw2+pW8lJi2bTUOl0dURE27ijP5rmadiad0j8+VffQOE09I2yw2CC4V35qLLERkjQCQZLGVH6ePeW1uVB2FLuR7qFxjlzqYt+GXFP3Y5jJvWty6B+d5OXLMgNuqlMt7vEMKw6Cg/vvuCDRJjXEAkCSxhT+LCMylffds4xrzOyps204XZp9663TNeW1qySDuCg7z56XLqqpqpr6sNvU9YV0VlSQaOPCtQFcLlkZ7k+SNKbwY5HbN8hMjKYwPY4TjbIyfCYHzrWxLCOelUsSzQ7lTWIi7ewuyeD5mg4plz5FdXMfpdmJxEVZd0uepUk2hiecNPWMzH2w8JkkjakCMHvKa/PSVF5r7JMXnmn6RiZ45XI3d69ZYrmuKa/bV2XT2j8mXR0eLpe1KtvOZmmi++VNxjX8S5LGFP6uPTXVpsJUuobGae4ZDcwFgtSz59txuDR3r7ZuSdJby7IAeP5Cu8mRWEN91zCDYw7Ljmd45SXYsClJGv4mSWOaQL3XrShyj2tIF9UbHTjXRl5KrCV2fZtNZmI06wtSeK5GSorAlMq2Fp055RVlVyzPTOC8tBD9SpLGFIFsaZRkJZIQHXH9D07A0LiDw3Vd3LXaul1TXntXZnGqpU/KwQBVTb0kRkdQnGmdyrazKc9JkpaGn0nSmMLf+2lMZbcp1uQlcUqSxnXPX+hgwuHi7jXW7Zryum1VFlrDoVppbVQ397GuIBlbEGxQtSoniat9o/SPTJodSsiQpDFFoKbceq0vSKHm2qAs8vM4cLaNjIRoNhdaZ0HfbMpzkshJjuH5MO+iGp1wcqFt0PKD4F6rctwz8mraQr+1MTg2icuAiTaSNKYLYDfJ+vwUJpwuLkgfK2OTTg7VdnDX6uyg2FJVKcVtK7N4sa4zrJP+2dZ+nC5tySKFMwmnciJ//nA1f/tK4HcJlaTh4Z0KG+iWBsDpFumieuFiJyMTzqDomvK6bWUWwxNOXg3jnRhfr2wbHC2NzMRo0uOjwiJpNPaMkBEb+DdgkjQ8vK26QI7H5ibHkJEQRXVzf+AuEiQOnG0jOTby+pa4wWDHinSi7DYOX+w0OxTTVDf3kZcSS2ZitNmh+EQpRdmSRGrbh8wOJaBcLk1TzwiZsYF/SZek4eHtCQzU4j5w/wKvz08J+5aGw+ni+doObluZRaRFyqD7Ii4qgoqiVA5fDN9S6VVNvZbbqW8uZUsSqWsfDOlyIm0DY0w4XGTFSUvDMNe7pwL8f74uP4VLnUMMjYfv5jBVzX30jUxy+6oss0OZt5tLM6ltH6StP/B9x1bTMTBGa/9Y0HRNeZVlJzIy4aSlN3QX1jZ0DwOQHSctDcO83tIIrPUFyWgNZ1rCt4vquZp2ImyKm0szzQ5l3m4uccd8uC78uqiqPNPFg62lUeqpaVbbHroTUJq63fW1pKVhIO+YRqDnnq/Ld//BnQrjLqqDNR1sW55GUkyk2aHM26qcRDITo8NyXKO6uY8Im2J1rnVX78+kNNudNC6GcNJo7Bkh0q5Ii5GkYRgj5jcDpMVHsTQtLmzHNRq7h7nUMcTtK7PNDmVBlFLcXJLJkUtdOEO4j3wm1U19rMpJIibSbnYo85IQHUF+amxIb/3a2D1MfmqcIdPXJWlMY0Q1i3X5yZwK0xlU3vpNwTie4XVzaQZ9I5OcuRo+z6HTpTndYv3KtrMpy07kYkgnjREK0+MMuZYkDY/rU24DPqrhnuN+tW80LOsYPX+hneKsBArT480OZcF2l2SiFGHVRVXXMcjwhDPoxjO8ypYkcrlziAlH6G3bq7V2J400SRqG0hgzewpeH9cIty6qgbFJjtX3BHUrA9xdjGvzksMqaQTbor7pypYk4nBprnQNmx2K3/UMTzA07mCpQW/EJGl4vN7SCLzVuUkoRVh1b4D7nbnDpdm7KjjHM6a6uSSTquY++kfDoxBedXMfybGRLMsIzhaidzA8FGdQNXp2JiyS7iljXZ9ya0DWiPeUlQ63abcHazpIjYu8vmd6MLu5NBOnS/PypfBY6Ffd3Mf6ghTLl7CfzYrMBCJsKiTHNRo9azRkTMNgr9eeMuaPYm1+Mqev9ofN9q8Op4tDtR3cWpYVFAUK57JxaQoJ0REcrgv9pDE07uBie/BUtp1JVISNZRnxITmDqrF7BKUgP1WShqGMbGkArMtLpnNwnPaB8BgMf30VePB3TQFE2m1sX57Gy5dDP2mcbunDpbH89q5zKV2SGJJrNRq7R8hJijFsKrQkDQ/tmVRhVPN7bZgNhntXge8uzTA7FL/ZWZxBY/cIzZ4+5VB1fXvXIE8aK7MTaeoZYWQitEr4NHYPGzobUZKGx/XZUwZdrzwnCbtNhc1geDCvAp/NrmJ3AnwpxMc1qpv6KEqPIzU+yuxQFsVbTuRiiFW8NXKNBpicNJRSdyulapVSl5RSn5/h+9FKqZ97vn9MKVUUqFiMKI0+VWyUnZKsBE6HwWB4sK8Cn01xVgJZidG8dLnb7FACRmtNdXPwLuqbqsxbTiSExjWGxh10D0+wNByShlLKDnwPuAcoB96tlCqfdtiHgV6tdTHwz8A3AhWPUQULp1qXn8yZMBgM964CD4WptlMppdhZnMHLl7pCtuz2tf4xOgbHQyJpLE2LIybSFlLTbr0zp4rCpHtqK3BJa12vtZ4AHgb2TztmP/Ajz+NfArerAA06vF4a3bi0sTY/hZ7hCa72hW7JZoCDNe2UZCUY+m7IKDuLM+gengipF6KpqryL+kJgmrTNpijNTqQ2hFoajZ7qtksNWg0OEGHYld4sD2ie8nkLsG22Y7TWDqVUP5AOvKETWSn1IPAgQHZ2NpWVlfMOZmDcnTTq6uqonGiY988vxGS/e6/pnx14mS1LzHkqhoaGFvT/5auRSc2x+hHuKooM6HXmw5/3bBtzz6D40dPHuHuZdcdrFnrPj18YJ0JBZ10VlZeDa6r0TPecpMc50+y0zO/iYh2qnwCgueYkXXUq4H/PYG7S8But9UPAQwAVFRV6z5498z6Hw+licKKSe2/bZdiA3w6Hk6+9egBXSj579qw05JrTVVZWspD/L189cboVp67iQ3dVsKUoLWDXmQ9/3/N3z1XSRhx79mz12zn9baH3/L0LL7O2QHPHbTv9H1SAzXTPl+z1HPldDeu23ERakA/sAzzdfZr0+Hbu2XsrEPi/ZzC3e+oqUDDl83zP12Y8RikVASQDARl1jLDbyEu0GTpDJDrCTtmSxJBeGR5Kq8Bns3NFBq9e6Qm5YniTThdnrvaHxHiG1/VyIiHSRWX0zCkwN2kcB0qUUsuUUlHAu4DHph3zGPABz+PfB57XITZqvDbPvWd4iN0WEHqrwGezsziDkQnn9fUMoaK2bZCxSVdIJY0yz7Tbuo7QSBpNPSOGV4w2LWlorR3AJ4ADQA3wC631OaXUl5VS+zyH/SeQrpS6BHwaeNO03GC3Lj+ZgTEHTSG4QOxkU2itAp/NjuXp2FTorde4vr1rQei0ErMSo0mOjQyJlsa4w0lr/6ihg+Bg8piG1vpJ4MlpX/vSlMdjwNuNjstIa/PcW2eebukP6j0mZnLwgncv8NBZBT6T5LhI1uYl89KlLj51R6nZ4fhNdVMfafFRFKTFmh2K3yil3BsyhcBst+aeUbSGoozw6Z4SuPtYoyJsnA3BleHeVeCJIbQKfDY7izOobu5jaDx0SlRUN/eyIYgr286mJDuB2rbBoO8S9q7RWJoWJt1Twi0qwsaqnKSQWxkeqqvAZ7OrOAOHS/PqldBYHd4/OsnlzuGgL1I4k7IliQyMOYK+WKh3jYZR+2h4+ZQ0lFLRvnxNLMy6vGTOXu0PqVXFoboKfDabClOJjrBxpC40koZ3UH9TYeiMZ3h5Z1AFexdVU88ICdERhk8d9rWl8YqPXxMLsDY/mcFxBw3dobMVZSivAp9JTKSdLUVpITMYXtXUi1LuiRqhJlSSRkP3MIXpcYZ3H94waSilliilNgOxSqmNSqlNno89QHi8GhjA+4cZKhVvB8YmefVKT8jPmppuZ3EGte2DdAyOmR3KolU19VGalRiS41Fp8VFkJkYH/QyqJhPWaMDcLY27gH/AvfDuH6d8fAr4q8CGFj6KMxOIibSFzLjGC7XevcCzzA7FUDuL0wF4Jcir3rpc7sq2G5eG3niGV7DPoHK6NM29xq/RgLmTRorW+lbgG1rr27TWt3o+9mutf21EgOEgwm5jdW5yyKwMP1jTTlp8FBtDeBX4TFbnJpMcGxn0XVRXuofpH50M6aRRkp3AxfahoB1HbO0bZdKpKTR4jQbMnTQ+5Pn3bYEOJNytzUvmbGs/ziD9JfZyrwLvZE9ZZkivAp+J3abYsTydI3VdQT2d01vZNpSTfll2IqOTTlp6g7PC9PXqthbsnqpRStUBZUqp01M+ziilThsRYLhYl5/MyIST+s7g3lXstcZe+kcnw2bW1HQ7SzJo7R/jSlfwTmqoauolMTqC4swEs0MJmNd38QvOLqrGHuP30fC64YpwrfW7lVJLcJf62HejY8XieAfDT7f0U+KZ3RGMDl7oINKu2F0S2qvAZ7N7yhawy4P0RbeqqY8NS1OwhXBLPxw4AAAeeklEQVRLsSTL/dzUtg+ytzz43uA0dY8QFWFjSVKM4deec8qt1rpNa71ea904/cOIAMPFsowE4qPsQT+D6rmadrYvTw/JWTe+KEyPIy8llhfrgnNcY2TCwYW2gZBc1DdVYkwkeSmxQdvSaOgeZmlanCmJ/YYtDaXUL7TW71BKneH1HVHBvSuq1lqvC2h0YcRuU6zOS+Z0S/BWSr3SNUx95zDv315odiimUcrdyvrdmWs4nC4i7MFVdOF0Sz8uHdrjGV6lnnIiwaixe8SUQXCYu6XxZ55/7wPeMuXD+7nwo3V5yZxrHcDhDM59GQ7WtAOE3fqM6XaVZDA45uB0ELYar2/vGuItDXCPa9R3DjMZZH9vWmvPPhrmFDi9YdLQWl/z/PumrinpnvK/tfnJjDtc1HUE52D4czXtlGUnUmDSOyCruGlFBkrBkSDsoqpq6mVZRryhm5GZpSw7kQmn63rhv2DRPjDO6KSTZQZXt/Waa0X4oFJqYIaPQaXUgFFBhot1+e53d8G4XqN/dJLjDb3cHmYL+maSFh/F6twkjgTZeg2tNSebQntR31SvlxMJrjdp3pl5RRnWbGkkaq2TZvhI1FonGRVkuChMiyMxJoLTV4NvXOOFi504XVqShseu4kyqmnoZDqJS6S29o3QNjYfFeAZAcVYCNhV8W796a9SZMd0WpDS6pdhsijVBujLcuwp8Qwjt8rYYu4ozmHRqjgVRqfSTTb0AIT9zyism0k5henzQzaBq6BomKsJGboo5m2NJ0rCYdfnJ1FwbZMIRPINzk04Xhy6E/l7g81FRFHyl0o839JAQHcGqnPDpRCjNTqA2yJJGfdcwhWlxpv2tSdKwmLX5yUw4XUH17udofTcDYw7uWh3es6amiom0s3VZGkcudZodis9ONPSyqTA1rBJ/WXYiDV3DjE06zQ7FZw1dw6aNZ4AkDctZl+fuGgimircHzrURG2nn5tJMs0OxlJ3FGVxsH6J9wPql0vtHJqltH2RLCG66dCOlSxJxaajvDI4ZVC6XprFnhGWSNIRXQVosybGRnAmSwXCXS/PMuXb2lGUSE2k3OxxL2TWlpIjVnWzqRWvYsizN7FAMVRZkGzK19o8y4XBJ0hCvU0qxLj85aFoaVc19dAyOc9fqJWaHYjnlOUmkxUcFxXqNVxt6iLQr1ueHxyC4V1FGPJF2FTTjGten25o0cwokaVjS2rxkatsGg6Kf9ZlzbUTYFLeulKm209lsiptWpHPkkvVLpZ9o6GFNXjKxUeHVWoy021iekcDFIJl22+BJGtLSEG+wLj8Fh0tzrtXarQ2tNQfOtXFTcQbJseFZoHAuu0sy6Bgct/Qq/7FJJ6ea+9lSFF5dU16lSxKDqKUxQmykneykaNNikKRhQRVF7sHIEw29JkdyY7XtgzR0j8isqRvYVeKeHGDlqrdnr/Yz4XRREWaD4F5l2Qm09I4GxULMhu5hCtPjUMq8GW6SNCwoIyGa5RnxHLd40jhwth2l4I4g3I/AKHkpsSzLiOdInXWn3np/zzaHadLwlhOxcmvQq6FrmOWZ5nVNgSQNy6ooSuW1xh5L72F84Fwbm5emkpVo/EYwwWR3SQZH63sYd1hzjOpEQw8rMuNJTzCvy8NMZd5d/Cw+ruFwumjqGTF1EBwkaVhWRVEavSOT1HdZ891Pc88I568NyKwpH9xSmsnopJPjV6zXcnS5NCcae8N2PAOgIDWOmEib5cc1WnpHcbi0qQv7QJKGZXn/iK3aRfXE6WsA3L1GksZcdqxIJyrCxqHaDrNDeZOatgH6RyfZtjx8k4bNpijJSrT8Wo0r3ebPnAJJGpZVlB5HRkIUxxt6zA5lRk+cbmVDQUrY753hi7ioCLYvT6fSgknjlcvu2lg7lofnnu5epdmJlq9222CBNRogScOylFJUFKZZcgZVfecQ51oHuG9djtmhBI09pZlc7hymqXvE7FDe4OXL3SzPjGdJcniPS5UtSaBjcJy+kQmzQ5nVla5hEqMjyEgwd4MsSRoWVlGUSlPPiOVqFz1x+hpKwX3rcs0OJWh4Fz9WXrROa8PhdPHqlR52LE83OxTTBcOGTFc8hQrNnG4LkjQszTuuceyKtbqoHj/VypbCtLB/dzofyzLiKUqPo7LWOlNvz1ztZ2jcwU0rwrtrCl5PGlYeDG/oNre6rZckDQtbk5dMYkwEL1uo4F1t2yB1HUO8Zb10Tc3XnrIsXr7cZZnyMK/Uu8cztofxILhXTnIMidERlp12Ozbp5GrvqOmD4CBJw9LsNsWO5emW2mv68VOt2BTcvUaSxnztKctkbNLF0XprbMz0yuVuVi5JDNv1GVMppSxdTqSxewSXhhUmL+wDk5KGUipNKfWsUqrO8++blqIqpTYopV5RSp1TSp1WSr3TjFjNtrM4g5beUUsMoGqteeJ0KztWpJOZKC8087V9eToxkTZLdFFNOFwcb+hhu4xnXFea7Z52a8Xikpc73WMtKzITTI7EvJbG54GDWusS4KDn8+lGgPdrrVcDdwPfVkqFV91m3EkD4KXL5rc2Tjb10tA9wv71eWaHEpRiIu3sWJ7OodoO01+Yqpv7GJt0cdMKSRpeZdkJ9I1M0jk0bnYob3LZU+LE7BIiYF7S2A/8yPP4R8Bbpx+gtb6ota7zPG4FOoCw2xpuRWY82UnRluii+uVrLcRG2rlXptou2G2rsmnsHrn+ztEsR+o6sSnYtkyShlfp9XIi1ptBdblziLyUWOKiIswOBbMiyNZaX/M8bgNuWPFOKbUViAIuz/L9B4EHAbKzs6msrFxQUENDQwv+2UAqTnByuOYazx86hM3P0+18vedxp+bRkyNsyorgxCtH/BqD0cx8nhPGXAD82+Ov8JYVxs23n37Pj58YZXmyjapXXzIsBqPN93keGHe3/n73UhWOq9Yq9V9dP0pqlJrzfgz53dZaB+QDeA44O8PHfqBv2rG9NzhPDlALbPflups3b9YLdejQoQX/bCD98kSzLvzcE/rc1X6/n9vXe360qkUXfu4J/dKlTr/HYDSzn+f93z2i9/2/Fw295tR77hwc04Wfe0J/57mLhsZgtIU8z5u+/Iz+y0dO+T+YRXC5XHrVXz+l/+axs3Meu5jfbeCE9uE1NmAtDa313tm+p5RqV0rlaK2vKaVycHc9zXRcEvA74Ita66MBCtXyvOMah+s6Kc9NMiWGX77WQl5KLNulO2PR7lydzTefruVa/yg5ybGGX/9FT5n2PWWy2+J0pdnWm0HVNjDGyITTEoPgYN6YxmPABzyPPwD8dvoBSqko4DfAj7XWvzQwNstZkhxDeU4Sz9eYs5q4tW+UI5e6eGBzPjabuatRQ8Gd5e4ij8+ebzfl+pW1nWQkRLHapDcgVla2JJE6i82gutzhrjkV7knj68AdSqk6YK/nc5RSFUqpH3iOeQdwM/BBpVS152ODOeGab++qLE409tA7bHxtnEdOtKA1PLBJZk35Q3FWAisy4zlwrs3waztdmsMXO7m5JFPeAMygNDuR4QknV/tGzQ7luuvTbbPMnzkFJiUNrXW31vp2rXWJ1nqv1rrH8/UTWuuPeB7/VGsdqbXeMOWj2ox4rWBveTYujeHltSedLv731UZuLs2k0OTqmqHkztVLOFrfQ//IpKHXPXO1n96RSW4pC7uJiD4pzXa/m7dSmfTLnUMkxkSQaZFFmLIiPEisyU0mKzGa52qM7dJ49nw77QPjvG97oaHXDXV3lmfjdGnDn89DFzpQCnaXSNKYSYm3BpWFpt1e7hxiRWaC6YUKvSRpBAmbTXH7qiwOX+xiwuEy7Lo/fqWBvJRYblspg6b+tD4/hbyUWB4/3WrodZ85387mpamkxZtbXtuqkmMjyUmOsVZLo2PYMuMZIEkjqOxdlc3QuMOw1eEX2wc5Wt/De7YvxS79335lsynuW5/Dkbouegwap2rqHqFGtuidk5U2ZBoad9A2MGaZ8QyQpBFUdpVkkBQTweOnjHl3+t8vNRBlt/HOigJDrhdu9q3PxeHSPHnm2twH+4F34F2Sxo2VLUnkUucQTpf5M6jqLVRzykuSRhCJjrBzz5ocDpxtC3h57Y6BMX71WgsPbM6XKqgBUp6TxIrMeB4z6E3AgXNtrMpJYmm6bNF7I2XZiUw4XFzpMn9c41KHN2lIS0Ms0L4NuQxPOHn+QmBnUf3XSw04XC7++OblAb1OOFNKsW99HscberjWH9gpnn3jLl5r6uWu1Tes2CPg+gLac60DJkfi3kkw0q4sNXNRkkaQ2b7cXZb80aqrAbvGwNgk/3O0kXvW5lhip7BQtm9DLloT8C7Hk+1OtJauKV8UZyUQZbdx3hJJY5AVmQlE2q3zUm2dSIRP7DbF2zbmcfBCR8D2Dv/vIw0Mjjv46C0rAnJ+8bplGfFsKEjxLKAMXB/6y60OSrMTWOmp5CpmF2m3UbokwRItjdq2wetb0VqFJI0g9O6tS3G6ND8/3uz3c3cPjfP9F+u5a3U2a/KS/X5+8Wbv2lJAXccQJ5v6AnL+pu4RLvW5uH9jvmXm+lvd6pxkzl8bMLWcyNC4g6t9o5RZLNFL0ghCRRnx7C7J4GevNuFw+nfNxr9WXmZkwsFn7izz63nF7O5bn0tclJ2fH28KyPl/4+nK3L8hNyDnD0Wr85LoGZ6gLUCteV/UedaKlGRZZ+YUSNIIWu/ZVsi1/jGe8WPRu+aeEX5ytJEHNuVfXxkrAi8hOoK3rMvl8VPXGBzzb1kRrTWPVl9lZZqN3BTjK+oGq/Ic92C4meMa3gWG0tIQfnFHeTbLMuL57vOX/NaE/tvHz2FXik/dUeqX8wnfvXNrAaOTTh6t9u+A+MmmXq50DbMj1/wd34LJypwklDJ3BlVt2xAxkTYKUq01RVqSRpCy2xQfv7WY89cGOOiHkunPnGvjuZoO/nxvibwjNcHGghTW5Sfz30eu4PLjorKfHm0iMTqCbUskacxHQnQERenxnGvtNy2Gi+3uQXCrVSOWpBHE9m/IpSAtlm8fvLio1asDY5P87ePnKctO5A93LfNjhMJXSik+vGsZ9V3Dfqtk3DU0zu9OX+OBzfnERFjrhScYlOcmcf6aiS2NduvNnAJJGkEt0m7jM3eWcfbqAL84sfCZVF969CxtA2N87YG1lpoPHm7uXZtDbnIMP3jxil/O94sTzUw4Xbx3+1K/nC/crM5NorlnlP5RY8vXA/QOT9A5OE6ZJA3hb/vW57K1KI1vPn2BjsH5z/R4sWWSR6tb+bPbS9i0NDUAEQpfRdptfHBnEa/Ud1PdvLjpt2OTTn74UgM7i9MpzrLeC08w8A6G15jQ2vAOgpdkW2vmFEjSCHpKKf7+bWsYnXTyF784Na/+8Fcud/PDcxPctCKdj99aHMAoha/+YFshafFR/OMztYs6zyMnmukYHOfje+R5XajVue51SmYMhlt15hRI0ggJxVmJfOm+1bxY18XXnqrxaTZVVVMvf/yTE2THKf7tvZul9LlFJERH8NFbVvBiXRdH67sXdI4Jh4t/f6GeTUtT2LEi3c8Rho/MxGiyEqNNGQyvaRskKSaCJUkxhl97LpI0QsS7txbwwZuK+P6LV/j6UxduODD+9Nk23vuDY6TERfHpihiSYyMNjFTM5X07CslOiuZrT9YsaILD/x5r5GrfKJ+8vURWgC9SeW4S564a39I43zpAeW6SJZ8/SRohQinFl+4r573bl/Ifh+v5g+8f5VRz3xtaHZc6hvizh6v4k5++xoqsBB75kx1kxMqvgNXERNr5q3tXcaqln/851jivn+0dnuCfn6tjV3EGe0plS9fFWpefQl3HICMTDsOu6XRpLrQNUJ5jzTI+Mnk7hNhsiq++dS3r8lL42lM17P/eS+SlxJKbEkPn4DgN3SPERNr45G3FfPK2EqIibNSYHbSY0b71uTxyooVvPl3L3lXZPq+d+ervahgad/DX95Vb8l1qsFmfn4xLw9mrA2xdlmbINa90DTE26WK1p0S71cjbzBD0ji0FVH72Vv7u/jVsWJqCTSnKc5P46/vKOfyXt/IXd5YRFSFPvZUppfjqW9egteaTP6ti0ocaY0+ducavTrbw8VuLLTmAGozW5acAcLolMMUkZ+IdeC+3aNKQlkaISo6N5D3bCnnPtkKzQxELVJQRz9cfWMcnf1bFF39zhq+/bd2sq4PPXu3nLx45xfr8ZD55m8yY8pfMxGjyUmIXPQV6Ps63DhBlt1FssUKFXpI0hLCwt6zPpa5jiO8crENr+Mpb1xATaX/DMa819vLgj0+QGhfFQ++vkAWafra+IJnTLcbNoDp/bYDSJdbaeGkqSRpCWNyn9pYA8J2DdZxu6efjtxWzuTCVvpEJHq26yg9fbiAvJZb//tBWsi04RTPYrctP4ckzbfQMT5AWHxXQa2mtOdc6wN5VWQG9zmJI0hDC4pRSfPqOUjYUJPO3j5/nT39WNeV7cP/GPP7698pJDfALWrhaP2VcY09ZYF/M2wfG6RmeuL4a3YokaQgRJG5bmc0tpVmcbOrlcscQcdERbC1KY0mytC4CaW1+MkrBqeb+gCcN70LC8lxrTrcFSRpCBBW7TbGlKI0tRcZM/xTuVfrFmQmGzKA63dKPTWHZ6bYgU26FEGJO6/JTqJ62WDYQqpv7KMlKJD7auu/nJWkIIcQcNhWm0D08QWP3SMCuobXmVEsfGwpSAnYNf5CkIYQQc6godHcHnmjsDdg1mnpG6BuZZL0kDSGECG4lWQkkxURwoqEnYNfwLiBcX2DdQXCQpCGEEHOy2RQVRWkBbWlUN/cRE2mz5G59U0nSEEIIH2wuTOVSxxA9wxMBOf+p5j7W5iUTYdGV4F7Wjk4IISzCO835tQC0NiYcLs62DlxfSGhlpiQNpVSaUupZpVSd599ZN6dWSiUppVqUUt81MkYhhJhqXX4ykXbFiUb/j2vUXBtgwuFiw1JJGrP5PHBQa10CHPR8PpuvAIcNiUoIIWYRE2lnbV4yr17xf9I47hlgD4ZFm2Yljf3AjzyPfwS8daaDlFKbgWzgGYPiEkKIWe1Ykc7pln4Gxyb9et5jV3ooSo8LioKTKtArHGe8qFJ9WusUz2MF9Ho/n3KMDXgeeC+wF6jQWn9ilvM9CDwIkJ2dvfnhhx9eUFxDQ0MkJFizhn2gyD2HB7ln/6jpdvKN42P8+aZoNmT5Z9W2S2s++fwIm7Ii+PDa6EWdazH3fOutt76mta6Y67iArVVXSj0HLJnhW1+c+onWWiulZspcHwOe1Fq3zLVtpdb6IeAhgIqKCr1nz54FxVxZWclCfzZYyT2HB7ln/9g+6eTbVc/QH5vDnj2r/XLO2rZBhg8cZv9Nq9mzOX9R5zLieQ5Y0tBa753te0qpdqVUjtb6mlIqB+iY4bAdwG6l1MeABCBKKTWktb7R+IcQQgRMTKSdrcvSeOlSl9/O+eqVbgC2GbQH+WKZNabxGPABz+MPAL+dfoDW+j1a66Va6yLgM8CPJWEIIcy2sziDi+1DdAyO+eV8R6/0sCQphvzUWL+cL9DMShpfB+5QStXhHq/4OoBSqkIp9QOTYhJCiDntKs4A4OVL3Ys+l9OleelSFzuLM5irG94qTEkaWuturfXtWusSrfVerXWP5+sntNYfmeH4H842CC6EEEYqz0kiNS6Sw3Wdiz7Xmav99I1McktZph8iM4asCBdCiHmw2RR7yrI4dKEDh9O1qHO9UNuJUrDb03oJBpI0hBBinu4oz6Z3ZHLRJUUO13WyLj8lqPZ3l6QhhBDzdHNpJlF2G8/VtC/4HP0jk1Q19XJLSfC0MkCShhBCzFtCdAQ7VqTz7Pn2BW8BW3mxA5eGW8qy/BxdYEnSEEKIBbijPJuG7hEutA0u6OefPHON7KRoNlp8p77pJGkIIcQC3Ls2hwib4tGqq/P+2eFxB5W1ndyzJgebLTim2npJ0hBCiAVIi49iT1kmj1ZfxemaXxfVodoOxh0u7lkzU6Ula5OkIYQQC3T/xnzaB8Z55fL8Fvr96rUWliTFUBEEpdCnk6QhhBALdPuqLBJjInjktWaff6atf4wXLnby+5vzsQdZ1xRI0hBCiAWLibTzwKZ8fnf6Gu0DvtWi+uVrzbg0vKOiIMDRBYYkDSGEWIQP7SzCqTU/fqVhzmPHJp386JVGdhVnsDQ9LuCxBYIkDSGEWITC9HjuLM/mJ6800jcyccNjf3Wyhc7BcT66Z4VB0fmfJA0hhFikP99byuC4g+8dujTrMWOTTv6t8jLr85O5aUW6gdH5lyQNIYRYpFU5Sfz+pnx+9HIjF9tnXuz37y9cpqV3lM/dvTJoyqDPRJKGEEL4wV/evZLEmAj+9GdVDI073vC9Ew09/L/nL7FvfS43BVFF25lI0hBCCD/ITIzmH9+xnrqOIT74X6/S2jcKwOGLnfzhD4+TnxrLV+9fY3KUixewPcKFECLc7CnL4l/etYHPPnKa3d88RHZiNK39YyzPjOfHf7iVpJhIs0NcNEkaQgjhR/ety2V9fgo/P95MS+8IG5em8s4tBcRE2s0OzS8kaQghhJ8VpMXxmbvKzA4jIGRMQwghhM8kaQghhPCZJA0hhBA+k6QhhBDCZ5I0hBBC+EyShhBCCJ9J0hBCCOEzSRpCCCF8prSe34boVqeU6gQaF/jjGUCXH8MJBnLP4UHuOTws5p4LtdaZcx0UckljMZRSJ7TWFWbHYSS55/Ag9xwejLhn6Z4SQgjhM0kaQgghfCZJ440eMjsAE8g9hwe55/AQ8HuWMQ0hhBA+k5aGEEIIn0nSEEII4bOwTBpKqbuVUrVKqUtKqc/P8P1opdTPPd8/ppQqMj5K//Lhnj+tlDqvlDqtlDqolCo0I05/muuepxz3gFJKK6WCfnqmL/eslHqH57k+p5T6X6Nj9DcffreXKqUOKaWqPL/f95oRp78opf5LKdWhlDo7y/eVUuo7nv+P00qpTX4NQGsdVh+AHbgMLAeigFNA+bRjPgb8u+fxu4Cfmx23Afd8KxDnefzRcLhnz3GJwGHgKFBhdtwGPM8lQBWQ6vk8y+y4Dbjnh4CPeh6XAw1mx73Ie74Z2AScneX79wJPAQrYDhzz5/XDsaWxFbikta7XWk8ADwP7px2zH/iR5/EvgduVUsrAGP1tznvWWh/SWo94Pj0K5Bsco7/58jwDfAX4BjBmZHAB4ss9/xHwPa11L4DWusPgGP3Nl3vWQJLncTLQamB8fqe1Pgz03OCQ/cCPtdtRIEUpleOv64dj0sgDmqd83uL52ozHaK0dQD+Qbkh0geHLPU/1YdzvVILZnPfsabYXaK1/Z2RgAeTL81wKlCqlXlJKHVVK3W1YdIHhyz3/DfBepVQL8CTwSWNCM818/97nJcJfJxKhQSn1XqACuMXsWAJJKWUD/gn4oMmhGC0CdxfVHtytycNKqbVa6z5TowqsdwM/1Fr/o1JqB/ATpdQarbXL7MCCUTi2NK4CBVM+z/d8bcZjlFIRuJu03YZEFxi+3DNKqb3AF4F9Wutxg2ILlLnuORFYA1QqpRpw9/0+FuSD4b48zy3AY1rrSa31FeAi7iQSrHy55w8DvwDQWr8CxOAu7BeqfPp7X6hwTBrHgRKl1DKlVBTuge7Hph3zGPABz+PfB57XnhGmIDXnPSulNgL/gTthBHs/N8xxz1rrfq11hta6SGtdhHscZ5/W+oQ54fqFL7/bj+JuZaCUysDdXVVvZJB+5ss9NwG3AyilVuFOGp2GRmmsx4D3e2ZRbQf6tdbX/HXysOue0lo7lFKfAA7gnnnxX1rrc0qpLwMntNaPAf+Juwl7CfeA07vMi3jxfLznbwEJwCOeMf8mrfU+04JeJB/vOaT4eM8HgDuVUucBJ/BZrXXQtqJ9vOe/AL6vlPoU7kHxDwbzm0Cl1M9wJ/4MzzjN/wUiAbTW/4573OZe4BIwAnzIr9cP4v87IYQQBgvH7ikhhBALJElDCCGEzyRpCCGE8JkkDSGEED6TpCGEEMJnkjSEWCSlVIpS6mOex7lKqV+aHZMQgSJTboVYJE/p/Ce01mtMDkWIgAu7xX1CBMDXgRVKqWqgDliltV6jlPog8FYgHnepjn/AXb77fcA4cK/WukcptQL4HpCJezHWH2mtLxh/G0LMTbqnhFi8zwOXtdYbgM9O+94a4G3AFuDvgBGt9UbgFeD9nmMeAj6ptd4MfAb4V0OiFmIBpKUhRGAd0loPAoNKqX7gcc/XzwDrlFIJwE28Xr4FINr4MIXwjSQNIQJrarVg15TPXbj//mxAn6eVIoTlSfeUEIs3iLvU+rxprQeAK0qpt8P1/Z3X+zM4IfxJkoYQi+SpEvuSUuos7mrB8/Ue4MNKqVPAOWbellYIS5Apt0IIIXwmLQ0hhBA+k6QhhBDCZ5I0hBBC+EyShhBCCJ9J0hBCCOEzSRpCCCF8JklDCCGEz/4/E4/qTwdyJTIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot drag/lift forces over time\n",
    "try:\n",
    "    import matplotlib\n",
    "    import numpy as np\n",
    "    import matplotlib.pyplot as plt\n",
    "    %matplotlib inline  \n",
    "\n",
    "    plt.plot(time_vals, drag_x_vals)\n",
    "    plt.xlabel('time')\n",
    "    plt.ylabel('drag')\n",
    "    plt.title('drag')\n",
    "    plt.grid(True)\n",
    "    plt.show()    \n",
    "\n",
    "    plt.plot(time_vals, drag_y_vals)\n",
    "    plt.xlabel('time')\n",
    "    plt.ylabel('lift')\n",
    "    plt.title('lift')\n",
    "    plt.grid(True)\n",
    "    plt.show()    \n",
    "except:\n",
    "    pass"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
